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Abstract 



The first attempts at solving a binary black hole spacetime date back to the 1960s, with 
the pioneering works of Hahn and Lindquist. In spite of all the computational advances 
and enormous efforts by several groups, the first stable, long-term evolution of the orbit and 
merger of two black holes was only accomplished over 40 years later, in 2005. Since then, the 
field of Numerical Relativity has matured, and been extensively used to explore and uncover 
a plethora of physical phenomena in various scenarios. 

In this thesis, we take this field to new frontiers by exploring its extensions to higher 
dimensions, non-asymptotically fiat spacetimes and Einstein-Maxwell theory. We start by 
reviewing the usual formalism and tools, including the "3+1" decomposition, initial data 
construction, the BSSN evolution scheme and standard wave extraction procedures. We then 
present a dimensional reduction procedure that allows one to use existing numerical codes 
(with minor adaptations) to evolve higher-dimensional systems with enough symmetry, and 
show corresponding results obtained for five-dimensional head-on collisions of black holes. 
Finally, we show evolutions of black holes in non-asymptotically flat spacetimes, and in 
Einstein-Maxwell theory. 
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Resumo 



As primeiras tentativas de evolugao da geometria de um sistema binario de buracos negros 
datam da decada de 60, com o trabalho pioneiro de Hahn e Lindquist. Apesar de todos os 
avancos computacionais c enormes esforgos por parte de varios grupos, as primeiras cvolucoes 
estaveis da orbita e coalescencia de dois buracos negros foram conseguidas apenas 40 anos 
depois, em 2005. Desde entao, o campo da Relatividade Numerica amadureceu, e tem sido 
usado extensivamente para explorar e descobrir fenomenos fisicos em varios cenarios. 

Nesta tese, levamos este campo a novas fronteiras e exploramos extensoes a dimensoes extra, 
espagos nao-assimptoticamente pianos e teoria de Einstein-Maxwell. Comegamos por rever 
o formalismo e ferramentas usuais, incluindo a decomposigao "3+1", construcao dc dados 
iniciais, o esquema dc evolugao BSSN c os procedimentos padrao de extraccao de radiagao. 
Scguidamente aprescntamos um proccdimento dc rcdugao dimensional que permite o uso de 
codigos numericos existentes (com adaptagoes menores) para evoluir sistemas em dimensoes 
mais elevadas com simetria suficiente, apresentando ainda os correspondentes resultados 
obtidos para colisoes frontais de buracos negros em cinco dimensoes. Finalmente, mostramos 
resultados em espagos nao-assimptoticamente pianos e teoria de Einstein-Maxwell. 
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Chapter 1 



Introduction 



Why you care about small things? World 
very simple place. World only have two 
things: things you can eat and things you 
no can eat. 

Quina 

Final Fantasy IX 



1.1 Motivation and historical background 

Formulated by Einstein in 1915 [HEIIS], general relativity is one of the most beautiful theories 
ever discovered. Its very elegance, however, can also be a disadvantage. We are able to do 
purely analytical calculations — sometimes using just pen and paper — in highly symmetrical 
ideal examples where exact solutions are known, or limits where gravity is weak. Alas, Nature 
is not that simple. 

To attack more complicated problems — such as those with strong and dynamical gravitational 
fields — one will eventually need to perform numerical computations. The quintessential 
example is the two-body problem. With well known solutions in terms of conies in Newtonian 
gravity, the equivalent problem in general relativity — the evolution of a black hole binary — 
posses no (known) closed-form solution. Perturbative analytical techniques do exist and some 
are very well suited to study certain stages of this problem. In particular, the inspiral phase 
(before the merger) is well modelled by post-Newtonian methods; the ringdown phase (after 
merger) can be described by perturbative methods using the quasi-normal modes of the final 
black hole. Full numerical simulations are required, however, to evolve the system during the 
merger. 

Much of the motivation to understand the nature of such systems and the corresponding 
energy emitted via gravitational radiation originally came from the gravitational wave as- 
tronomy field. A first generation of highly sensitive gravitational waves detectors — LIGO 
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Virgo [S] , GEO ^ and TAMA [7J — have been operational and a second generation of even 
more sensitive detectors is under construction. These detectors may allow us to study signals 
produced from strong-field systems, which carry the specific signature of the system that 
produced them. The analysis of these signals may then provide us with a new window to the 
universe. For that, however, we rely on source modelling: templates of theoretical waveforms 
from likely sources are needed if one wishes to reconstruct the signal. 

Numerical relativity can be regarded as a tool to study spacetimes that cannot be studied 
by analytical means. It dates back to the mid 1960s with Hahn and Lindquist's attempts of 
numerically evolving Einstein's field equations for a binary black hole spacetime [8j. Their 
computer power was very limited, however, and not much physics could be extracted from the 
simulation. More reliable simulations were only performed in the late 1970s by Smarr ^ and 
Eppley [To], which again attempted the head-on collision of two black holes. Though almost 
a decade after Hahn and Lindquist, the available computer power was still only enough to 
evolve low resolution simulations. 

With the development of faster computers and the extra motivation of returning to the two- 
body problem coming from LIGO, the 1990s finally witnessed the simulation of a head-on 
collision of two black holes [UllIS] as well as the study of more complex systems. To name just 
a few: simulations of rapidly rotating neutron stars [13], the formation of a toroidal event 
horizon in the collapse of a system containing a toroidal distribution of particles |14| [T5] 
and one of the most influential results, gravitational collapse and its relation with critical 
phenomena [16j. For a more comprehensive overview see, for example, |17| . 

In spite of all these successes, the real breakthrough came only in 2005 with the first simula- 
tions of stable, long-term evolutions of the inspiral and merger of two black holes |18 [ ll9 t l20j. 
For an overview of the two-body problem in general relativity see |21j . 

Since then, numerical codes have considerably improved and much progress has been made. 
In particular, we have witnessed numerical evolutions of (see e.g. [22] for a thorough overview 
of some recent results): 

• binary black hole mergers lasting for 15 orbits before merger |23ll24j . The corresponding 
waveforms, which include the infall, non-linear merger and ringdown phase, have been 
used in comparisons with post-Newtonian results. 

• black hole binaries with mass ratios up to 1 : 100 \25\ 126] . Waveforms for high mass 
ratios are essential, since they are expected to be the most common, astrophysically. 
Computationally, however, they are much more demanding than comparable mass cases. 

• rotating black holes with near extremal spin |27l [28] . 

• zoom-whirl orbits — characterised by black hole trajectories that alternate between a 
whirling quasi-circular motion and a highly eccentric quasi-elliptical zooming out [29^ 

ED]. 

• so called "superkick" configurations — equal mass black holes with (initially) opposite 
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spin vectors lying in the orbital plane — where the post-merger recoil velocity can reach 
up to ~ 4000 km/s |31i|32|. Recently [33], this effect was combined with the "hangup" 
configuration |34j — where the black holes have spin aligned with the orbital angular 
momentum — to predict maximum recoils up to ~ 5000 km/s. Given that such velocities 
are enough to eject a black hole from the centre of a galaxy [35j, these results could have 
important consequences for astronomy (such as in structure formation), see e.g. [MlEZl 

MM- 

• high velocity collisions of black holes — head-on collisions up to 0.94c [IQ], where the 
radiated energy was found to be around 14% when extrapolating the relative velocity 
between the black holes to c; as well as non-head-on collisions [411 US] ) where the impact 
parameter for black hole merger was determined in the limit where the relative velocity 
approaches c. 

As we can observe from these previous examples, numerical relativity has now reached a state 
of maturity and, at least in the four-dimensional asymptotically flat vacuum case, is largely 
under control allowing us to evolve a large class of different configurations. 

Motivation to study gravity in the dynamical/strong field regime is not restricted to the 
evolution of the two-body problem or variations thereof, and computation of the respective 
waveform. Indeed, incentive to study such systems also comes from fields other than gravity 
itself. In the following we mention some of these topics and briefly describe how numerical 
relativity can be expected to shed light on some outstanding issues. 

Tests of cosmic censorship hypothesis 

It has been known for some time from the Penrose-Hawking singularity theorems that, 
quite generically, solutions of Einstein's field equations with physically reasonable mat- 
ter content can develop singularities |l3]. If such singularities are visible to the rest of 
spacetime (i.e., no horizon is covering them), predictability may break down. Originally 
formulated by Penrose in 1969 [H], what is known as the weak cosmic censorship 
conjecture roughly states that, generically, singularities of gravitational collapse are 
covered by an event horizon and therefore have no causal contact with distant observers. 

In the absence of a generic prooQ one can try and put the conjecture to the test 
in specific configurations. The ability to perform full blown non-linear numerical 
simulations in arbitrary spacetimes could here prove invaluable. 

With such simulations, the conjecture has been shown to hold under extremely violent 
events — the ultra-relativistic collision of black holes [IQ] . In higher-dimensional gravity, 
it was shown by Lehner and Pretorius that cosmic censorship does not seem to hold 
generically, even in vacuum |45) . Specifically, it was shown that five dimensional black 
strings (solutions of five dimensional vacuum gravity, known to be unstable ) display, 
when perturbed, a self-similar behaviour that ultimately gives rise to naked singularities 

*Indeed, the conjecture has not even been stated in a rigorous way — as often happens, that is part of the 
task. 



CHAPTER 1. INTRODUCTION 



17 



in rather generic conditions. Also in higher dimensions, results by Okawa et al. [37] seem 
to indicate that in a high- velocity scattering of five-dimensional black holes, curvature 
radius shorter than the Planck length can be observed (i.e., no horizon is covering this 
region). This could be regarded as an effective singularity in classical gravity. 

Stability of (higher-dimensional) black hole solutions 

Higher dimensional gravity has a much richer diversity of black object solutions than 
its four dimensional counterpart. Spherical topology is not the only allowed topology 
for objects with a horizon and one can also have, e.g., black rings (with a donut-like 
topology) [IS HH] and even regular solutions with disconnected horizons, such as the 
"black Saturn" [50], the "black di-ring" [5T] or the "bicycling black rings" [52]. 

The study of these objects is relevant for a number of reasons. Other than the obvious 
intrinsic value that such studies carry and the possibly interesting mathematical prop- 
erties that some of these objects may have, the understanding of these solutions can 
also be helpful for: (i) quantum gravity — the calculations of black hole entropy within 
string theory were first performed in five dimensional spacetimes, and only afterwards 
extended to four dimensions; (ii) gauge/gravity correspondence, which maps properties 
of D-dimensional black holes to strongly coupled field theories in D — 1 dimensions; 
(iii) large extra dimensions scenarios, suggesting that (microscopic) higher-dimensional 
black objects could be formed in particle collisions with centre of mass energy > TeV 
(such as at the LHC). See the review article [53] for more details and further motivation. 

The stability of these higher-dimensional black objects is now starting to be explored. 
It had been conjectured that for D > 6 ultra-spinning Myers-Perry black holes will 
be unstable [SJ], and this instability has been confirmed by an analysis of linearised 
axi-symmetric perturbations in D = 7, 8, 9 [55j. 

Clearly, the study of the non-linear development of these instabilities and determination 
of the respective endpoint requires numerical methods. Such studies have been recently 
presented for a non axi-symmetric perturbation in D = 5 [56] and D = 6,7,8 j57] . 
where it was found that the single spinning Myers-Perry black hole is unstable, for 
sufficiently large rotation parameter. 

AdS/CFT correspondence 

In 1997-98, a powerful tool known as the AdS/CFT correspondence (or the gauge/gra- 
vity duality, or even, more generically, as holography) was introduced |58j- This holo- 
graphic correspondence (if true in general) is extremely powerful since it maps the 
dynamics of a non-perturbative, strongly coupled regime of certain gauge theories in D 
dimensions to {D + l)-dimensional classical gravity. This means that, for such gauge 
theories, we can map strongly coupled quantum field theory dynamics to systems of 
partial differential equations, which can in principle be solved (numerically, if needed). 

In particular, high energy collisions of black holes are said to have a dual description in 
terms of high energy collisions with balls of de-confined plasma surrounded by a confin- 
ing phase. These are the type of events that may have direct observational consequences 
for the experiments at Brookhaven's Relativistic Heavy Ion Collider (RHIC) \59\ 160] . 
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Numerical relativity in anti-de Sitter (AdS) is notoriously difficult, and therefore pro- 
gress has been slow in applying its techniques to the aforementioned problems. Never- 
theless, some exciting results have been recently put forth [6 H l62 l l63]. 

TeV-scale gravity scenarios 

As first pointed out by 't Hooft [64J , if two point particles collide at energies above the 
Planck energy, it is expected that gravity should dominate the interaction and thus, 
quite remarkably, the process should be well described by general relativity. 

Thome's Hoop Conjecture [65] further tells us that if one traps a given amount of 
Energy E in a region of space such that a circular hoop with radius R encloses this 
matter in all directions, a black hole is formed if its Schwarzschild radius Rs = 



recently gained more support with the work by Choptuik and Pretorius [66], where it 
was shown that collisions of boson stars do form black holes, for sufficiently high boost 
parameter. 

If this conjecture does hold, it would imply that particle collisions could produce black 
holes \67\ [68] • As argued above, the production of black holes at trans-Planckian 
collision energies (compared to the fundamental Planck scale) should be well described 
by using classical general relativity (see also |69] and references therein). Writing down 
the exact solution describing the collision of two ultra-relativistic particles in general 
relativity is not feasible, however, and approximations have to be used. One possible 
approximation (good for its simplicity) is to use black holes, and model the scattering 
of point particles by black hole collisions. 

This gains further relevance in the context of the so-called TeV-gravity scenarios. Such 
models were proposed as a possible solution to the hierarchy problem, i.e., the relative 
weakness of gravity by about 40 orders of magnitude compared to the other fundamental 
interactions. It has been proposed that this can be resolved if one adopts the idea that 
the Standard Model is confined to a brane in a higher dimensional space, such that the 
extra dimensions are much larger than the four dimensional Planck scale (they may be 
large up to a sub-millimetre scale) \70\ f7T\ [72] . In a different version of the model, the 
extra dimensions are infinite, but the metric has an exponential factor introducing a 
finite length scale [731173] . 

In such models, the fundamental Planck scale could be as low as 1 TeV. Thus, high 
energy colliders, such as the Large Hadron Collider (LHC), may directly probe strongly 
coupled gravitational physics [TSj jTB] [68l [671 [13 EHj- In fact, such tests may even be 
routinely available in the collisions of ultra-high energy cosmic rays with the Earth's 
atmosphere \79\ [80l I81| . or in astrophysical black hole environments \82\ l83l (for 
reviews see [HSl [13 [HH] ) • 

Numerical simulations of high energy black hole collisions in higher dimensional space- 
times, then, could give an accurate estimate of the fractions of the collision energy 
and angular momentum that are lost in the higher-dimensional space by emission of 




> R. This conjecture (or rather, the classical variant thereof) has 
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gravitational waves; such information would be extremely important to improve the 
modelling of microscopic black hole production, and of the ensuing evaporation phase, 
which might be observed during LHC collisions. 

The challenge is then to use the classical framework to determine the cross section for 
production and, for each initial setup, the fractions of the collision energy and angular 
momentum that are lost in the higher dimensional space by emission of gravitational 
waves. This information will be of paramount importance to improve the modelling of 
microscopic black hole production in event generators such as Truenoir [HZ], Charyb- 
Dis2 [87J, Catfish [88j or Blackmax [891 [90]. The event generators wih then provide 
a description of the corresponding evaporation phase, which might be observed during 
LHC collisions. 

For a thorough review of these topics, challenges and how tools coming from numerical 
relativity can help see |91j . 

1.2 The new frontiers 

With these motivations in mind, we propose in this thesis to extend current numerical 
relativity tools to new frontiers. 

1.2.1 Higher-dimensional gravity 

The first such frontier, in light of our discussion in the previous section, is higher-dimensional 
gravity. We start by emphasising that full blown 4 + 1, 5 + 1, etc. numerical simulations 
of Einstein's field equations without symmetry are currently (and in the near future) not 
possible due to the heavy computational costs. We have thus developed a framework and 
a numerical code that can, in principle, be applied to different spacetime dimensions (with 
enough symmetry) with little adaptations. This is achieved by taking the D dimensional 
vacuum spacetime to have an isometry group fit to include a large class of interesting 
problems. If this isometry group is sufficiently large, it allows a dimensional reduction 
of the problem to 3+1 dimensions, where our originally higher-dimensional problem now 
appears as (four dimensional) general relativity coupled to source terms. Thus, the different 
spacetime dimensions manifest themselves only in the different "matter" content of the four 
dimensional theory. An obvious advantage of this approach is that we can use existing 
numerical codes with small adaptations: taking a working four-dimensional code, the four 
dimensional equations need to be coupled to the appropriate source terms and some issues 
related to the chosen coordinates must be addressed. Incidentally, other issues possibly 
related with the choice of gauge conditions further complicate the problem. 

We should here point out that alternative approaches have been proposed to evolve numeri- 
cally Einstein's equations in higher dimensions, as well as other codes tailored to study specific 
problems. In particular we note the previously mentioned pioneering works concerned with 
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the non-linear development of the Gregory-Laflamme instability [46j of cosmic strings |92tl45j: 
studies of gravitational collapse [MlEl!; static situations [95j; and the alternative approach, 
based on the Cartoon method [96], that has been developed and tested by Yoshino and 
Shibata [97^ [57] . For a review of numerical relativity in higher-dimensions see also |98[ l99] . 

1.2.2 Non-asymptotically flat geometries 

Another frontier has to do with numerical evolutions in non-asymptotically flat spacetimes. 
de Sitter 

Going back to four dimensions, an obvious first choice of a non-asymptotically flat 
spacetime is de Sitter, the simplest model for an accelerating universe, de Sitter is 
a maximally symmetric solution of Einstein's equations with a positive cosmological 
constant, describing a Friedmann-Robertson- Walker (FRW) cosmology with a constant 
Hubble parameter. There is now a large body of observational evidence for a present 
cosmological acceleration well modelled by a positive cosmological constant A [100] . 

Cosmological dynamics should leave imprints in gravitational phenomena, such as 
gravitational radiation emitted in a black hole binary coalescence. Identifying such 
signatures can thus be phenomenologically relevant in view of ongoing efforts to directly 
detect gravitational radiation. 

Studying the dynamics of black holes in asymptotically de Sitter spacetimes can also 
potentially teach us about more fundamental questions such as cosmic censorship, as in 
the following scenario. Consider two black holes of sufficiently large mass in a de Sitter 
spacetime. If, upon merger, the final black hole is too large to fit in its cosmological 
horizon the end state of such an evolution would be a naked singularity. This possibility 
begs for a time evolution of such a configuration, which we will show and discuss. 

Black holes in a box 

As argued above, having a framework to solve Einstein's equations in asymptotically 
Anti-de Sitter geometries would be of major help for studies of AdS/CFT duality, in 
particular in dynamical settings. This is no easy task, however, and a major reason for 
that is the "active role" played by the boundary of AdS spaces. This is easily visualised 
in the Penrose diagram of AdS, which has a timelike boundary. Null geodesies in AdS 
reach the boundary in a finite affine parameter, and one therefore often refers to an 
asymptotically AdS space as a "box", having in mind that AdS boundary conditions 
directly affect the bulk physics [Ml [Ml EOH] . 

As a first step to model the role of the boundary in evolutions, we will here give an 
overview of the work reported in [104] where a toy model for AdS was considered. 
Therein, as we will explain, the cosmological constant is set to zero and mirror-like 
boundary conditions are imposed on a box containing the dynamical system, which 
mimics the global AdS geometry. 
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Inside this box, black hole binaries are evolved, producing an inspiralling merger. Such 
systems are very well tested in standard asymptotically flat spacetimes with purely 
outgoing boundary conditions, and differences to these cases can be clearly seen. 

Black holes in cylinders 

Again in the topic of higher-dimensional gravity, now in scenarios with compact extra 
dimensions, a natural question to ask is how the compactness of the extra dimensions 
changes the dynamics of such scenarios (as opposed to the asymptotically flat cases) and 
understanding the role of the compactness of the extra dimensions in the aforementioned 
TeV gravity scenarios. 

There is considerable literature on Kaluza-Klein black holes and black holes on cylin- 
ders \W5\ 11061 1107| [108] ■ but, to the best of our knowledge, the full non-linear dynamics 
of black holes in such spacetimes remain unexplored. 

In this spirit, using the formalism developed for higher-dimensional spacetimes, we have 
started exploring what happens when one of the directions is compactified. 

1.2.3 Einstein-Maxwell 

Finally, we have started exploring the electrovacuum Einstein-Maxwell system. We first note 
that while the dynamics of black holes interacting with electromagnetic fields and plasmas 
have been the subject of a number of numerical studies (e.g. jlOQllllO] ). dynamics of charged 
(Reissner-Nordstrom) black holes have remained rather unexplored. 

Studying the dynamics of charged black holes is relevant for a number of reasons. In the 
context of astrophysics, charged black holes may actually be of interest in realistic systems. 
First, a rotating black hole in an external magnetic field will accrete charged particles up to a 
given value, Q = 2BqJ |lllj . It is thus conceivable that astrophysical black holes could have 
some (albeit rather small) amount of electrical charge. It is then of interest to understand the 
role of this charge in the Blandford-Znajek mechanism jll2], which has been suggested for 
extracting spin energy from the hole, or in a related mechanism capable of extracting energy 
from a moving black hole [llOt I113j to power outflows from accretion disk- fed black holes. 
Also of interest is investigating the role of charge in post-merger recoil velocities of black hole 
binaries, and see if the recently predicted recoils of ~ 5000 km/s [33] could be exceeded. 

Incentive to study such systems also comes from outside of astrophysics. In particular, as 
already mentioned above, it was argued by 't Hooft [M], that in trans-Planckian particle 
collisions, gravity should dominate the interaction and thus the process should be well 
described by general relativity — we can say that, for ultra high energy collisions, matter does 
not matter [66j. Calculations of shock wave collisions, however, seem to suggest that even 
though other interactions — say charge — may become irrelevant in the ultra-relativistic limit, 
the properties of the final black hole (and of the associated emission of gravitational radiation) 
will in fact depend on the amount of charge carried by the colliding particles |1141 1115j . One 
then wonders whether the often repeated matter does not matter scenario is actually true. 
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Light can be shed in this issue by performing highly boosted coUisions of charged black 
holes (analogous to the ones performed in vacuum I41j ) and comparing the results — in 
particular the profile of the corresponding waveform — against equivalent electrically neutral 
systems. 

With these incentives, we will report on the first steps taken in the numerical evolution of 
Reissner-Nordstrom black holes, building on previous numerical evolutions of the Einstein- 
Maxwell system [HHl EMI IHZl EE] • 

1.3 Structure 

The structure of this thesis is as follows. We start by reviewing, in chapter [2j standard 
differential geometry results, summarise the formalism of the "(D— 1)+1" decompositiorj^Jand 
conformal decomposition to write Einstein's equations in a dynamical systems form. In chap- 
ter|3]we then review the construction of relevant initial data for the class of problems we will be 
interested on and discuss, in chapter[4j the numerical implementation of Einstein's equations: 
first, we need to re- write the evolution equations in the so-called BSSN (Baumgarte, Shapiro, 
Shibata, Nakamura) form; we next discuss the gauge conditions and finish by giving a very 
brief overview of the numerical code we use for the simulations. In chapter [5] we review the 
usual procedures to extract the physical results from numerical simulations: wave extracting 
and horizon finding. These chapters consist mostly of review material found in the usual 
literature (e.g. \119\ 1120] I121| ). Finally, in chapter [6| we introduce a dimensional reduction 
procedure that allows us to reduce higher-dimensional systems (with enough symmetry) to 
effective four-dimensional theories with source terms. This enables us to perform numerical 
evolutions of such higher-dimensional systems by adapting existing numerical codes. We 
also discuss the construction of initial data and present results. Chapter [7] is dedicated to 
evolutions in non- asymptotically flat spacetimes: we present the aforementioned collisions of 
black holes in asymptotically de Sitter spacetimes, black holes in a box and black holes in 
asymptotically cylindrical spacetimes. In chapter [8] we report on evolutions of charged black 
holes, in electrovacuum Einstein-Maxwell theory, and we end with some final remarks and 
future directions in chapter [9} 

1.4 Preliminaries 

Let us consider a D-dimensional pseudo-Riemannian manifold {M,g), that is, a differentiable 
manifold equipped with a smooth, symmetric metric tensor g with signature ( — !-•••+). 
We further assume that the manifold is covered by a set of coordinates {x^}, ^ = 0, . . . , D — 1. 

A coordinate basis of the tangent space of A4 at p, TpAi, is given by 5^ = d/dx^. A vector 

* Usually found in the literature as "3+1" decomposition. Here we will keep the D arbitrary, but the 
differences are minimal. 
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V € TpAi can be written in the form 

V = V^'di,, (1.4.1) 

where F'* are the components of V in the basis 9^. When a vector V acts on a function / it 
produces the directional derivative of / along V: 

Vif) = V''d^f. (1.4.2) 

A 1-form oj G T^M (the cotangent space at p) is an object which is dual to a vector, i.e., it 
produces a number when acting on a vector. The simplest example of such an object is the 
differential d/ of a function /. The action of d/ on V is defined to be 

{df,V) = V{f) = V^^d^f. (1.4.3) 

Since d/ = 9^/dx^, {dx^} is a natural choice as a basis of T^M.. We thus naturally expand 
an arbitrary 1-form oj as 

oj = uj^dxi'. (1.4.4) 

The metric tensor g allows us to define a scalar product between two vectors U and V 

U-V = g{U,V)=g {d^, 8,) C/^F^ = ^^.C/'^F'^ , (1.4.5) 

which induces an isomorphism between vectors and 1-forms, corresponding in the index 
notation to the usual raising and lowering of indices: if C/ is a vector, one can define a 1-form 
U\, through 

{U^,V) = g{U,V)=g^,U>'V'' = {U^)^V'' VF. (1.4.6) 
Analogously, given a 1-form oj we can map it to a vector a;" through 

{a, J) = g-\a,co) = g'^ (a^dx^, c^.dx'^) = g'^Vo;. = (uj^Y (l-^-^) 

Since wc will be working mostly in the index notation, and the placement of the index makes 
clear whether one is dealing with vectors of 1-forms, we will omit the flat and sharp symbols 
throughout. 

The metric further allows us to determine the distance between two nearby points in the 
manifold according to 

ds^ = g^^dx^dx" . (1.4.8) 

Notice that a basis of TpAi (and of TpAi) need not be coordinate. One can have, for instance, 
the combination = Aa^dfj,. {cq} is an example of a non-coordinate basis. 

We now introduce a (generic) non-coordinate basis obeying 

[ea,e/3] = Cc,/e5, (1.4.9) 
where the Lie bracket [X, Y] is defined by 

[X, Y]f = X {¥{/)) - Y {X{f)) . (1.4.10) 
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The connection coefficients T"^^ then take the form 

r"^7 = -g'^^ {gSPa + 557,/3 - 5/37,5 + C5/37 + C57/3 - Cp^i) , (1.4.11) 

where r"[^^] = — ^c^^". When using coordinate basis (0^7° = 0), these are usuahy called 
the Christoffel symbols. 

We now define the Riemann curvature tensor 

-^",375 = r"/35,7 - T^p^^s + r"A7r^/36 - ^"xS^^/B-y - /3\C^5^ ■ (1.4.12) 

Mind the notation 

Ta,p = de.Ta = efi (r„) . (1.4.13) 

Thus, Ta,aX = de.de^Ta ^ de^de^T^. 

General relativity is a geometric theory of gravity which relates the curvature of spacetime 
to its matter content via the Einstein field equations, which read 

Gf,, = Rf,, - ]^Rg^^ = SttT^, , (1.4.14) 

where = R^^xu is the Ricci curvature tensor, R its trace (the Ricci scalar), g^,^ the metric 
tensor and T^i, the stress-energy tensor. All these quantities are I?-dimensional. 

Throughout this work we will always use the ( — h h) metric signature and geometrised 

units G = 1 = c. 



Chapter 2 

{D — 1) + 1 decomposition 



We start by briefly stating some known results from differential geometry that will be of use 
to us. In this chapter we use the following conventions: Greek indices (a, /?, 7, . . . ) run from 
to D — 1; Latin indices {i, j, k, . . .) run from 1 to D — 1. 

We work on a D-dimensional manifold with a metric gf^^. We denote the torsion- free 
Levi-Civita connection on A4 associated with (7^,^ by ^ . All quantities defined relative to 
the manifold Ai will have a leading ^ superscript, e.g., the Riemann curvature tensor on Ai 
is denoted by ^R'^q/^^. 



2.1 Hypersurfaces 
2.1.1 Definition 

A codimension 1 hypersurface S is a (D — l)-dimensional submanifold of A4, defined as the 
image of a (-D — l)-dimensional manifold S by an embedding S = |122| 1119] . Given 

a scalar field t on A^, we can select a particular hypersurface S by putting a restriction on 
the coordinates 

t(x") = 0, (2.1.1) 

or by giving parametric equations 

= (yi) ^ (2.1.2) 
where are coordinates intrinsic to the hypersurface. 
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2.1.2 Normal vector 

The 1-form d^t is normal to the hypersurface. We can introduce a unit normal (if the 
hypersurface is not null) as 

ria = , ^ dgt. (2.1.3) 



With this definition 



{— 1 if S is spacelike 
(2.1.4) 
+1 if S is timelike 



2.1.3 Induced metric 

We obtain the induced metric on S by restricting the line element to displacements confined 
to the hypersurface. Using the parametric equations = (y*) , we define the vectors 

that are tangent to the curves in S. For displacements confined to S we have 

ds| = gi^udx'^dx" 

= jijdy'dy^, (2.1.6) 

where 

is the induced metric of the hypersurface (also called the first fundamental form ofE). Notice 
that if u, -y G S, 

u-v = gfxvU^v"" = ^iju^v^. 



2.1.4 Orthogonal projector 

The orthogonal projector onto S is a concept closely related with that of the induced metric. 
For a hypersurface E with unit normal n'* we define it as 

Pixu = gnu - crn^n^. (2.1.8) 

To see that this definition makes sense, we note that, for any vector in M. (or, more cor- 
rectly, in TpM., the tangent space of M. at p), P^,/ will project it tangent to the hypersurface, 
i.e., orthogonal to n^^: 



(P^.i;'^)n^ = 0. 



(2.1.9) 
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Notice also that the projector is idempotent: 

Pf'xP^^ = P\. (2.1.10) 

Finally, note that for n, f € E, P^y acts just like the metric: 

Pt^uu^'v'' = g^.uu^v" = j.ju'vK (2.1.11) 

Thus, we see that the orthogonal projector P^j^ is the natural extension of the induced metric 
7jj for all vectors in TpAi. As such, from now on we will no longer make any distinction 
between these two concepts, and will denote both by 'y^i, (defined as 7^,^ = g^u — crn^n^). 
This way we adopt a D-dimensional point of view, where we treat all tensor fields defined 
on E as if they were defined on Ai and we avoid the need to introduce a specific coordinate 
system on E. 

Note that we can project an arbitrary tensor on A4 onto E in the following way. Let 
T^^^"'^^p y^...y^ be a tensor field on M.. Denoting {'^T)°^'^"'"'^ another tensor in M. such 

that 

(7r)"--"- = 7"V, • • • ■ ■ ■ 7'''/3,r^^-^^.i....„ (2.1.12) 

we easily see that (7T)"^"'"^ I3i-i3q is in E. 



2.1.5 Intrinsic curvature 

We now want to define a covariant derivative associated with 7^1^ on E, V, that has the 
"usual properties" of a covariant derivative, in particular that it is torsion- free and satisfies 

Va7^^ = 0. (2.1.13) 



The easiest way to define it is just to project the covariant derivative ^ onto E using (2.1.8): 

Vp-t Pl--0q — I fll I tJ.pl I3l 1 Pql p V(j-t VI— Vq- ^^Z.l.i'lJ 

It can be shown jll9j that this definition of the covariant derivative has all the properties we 



want (linearity, Leibniz' rule, its torsion vanishes, . . . ) and it satisfies (2.1.13). 

We can now define the Riemann tensor associated with this connection, i?"^^^, as the measure 
of the non-commutativity of this covariant derivative, associated with the metric on E, 

V«V/3r;^ - VpVaV^ = R^apyf', (2.1.15) 

where £ T,. 

p^yS represents the intrinsic curvature ofT,. 
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2.1.6 Extrinsic curvature 

The intrinsic curvature of the hypersurface S, as the name impHes, is a property of the 
hypersurface itself. We wiU now define the extrinsic curvature, which depends on how S is 
embedded on M. We define it a^ |123) 

Kf,u = -l''^.l^yyanp. (2.1.16) 

It can be shown that K^^, = Kyfj_. Defining 

a^' = n^Vun^", (2.1.17) 

we have, after some simple algebra, 

Kfiu = -"^firiu + CFn^Qy. (2.1.18) 

We will always consider spacelike hypersurfaces, so from now on we will work with a = —1. 
Note that, by definition, K^y lives on E. 

2.2 Foliations 

We assume that our spacetime can be foliated by a family of spacelike hypersurfaces , that 
is, there exists a smooth scalar field t on such that 

^t = {p^ M, tip) = t} . (2.2.1) 

In the following we will not distinguish between t and t. 



2.2.1 The lapse function 

We write the unit normal vector as 



where 



= -ad^t, (2.2.2) 



a = , I ^ (2.2.3) 



is called the lapse function. 

*Our definition, with the minus sign, agrees with the standard convention used in the numerical relativity 
community, but note that some authors use different conventions. 
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2.2.2 Normal evolution vector 

We define the normal evolution vector as 

m^ = an^. (2.2.4) 

We can easily see that 

mi'^^t = 1, 

which means that m^, unhke n^, is adapted to the scalar field t. It can be shown |119j that 
the hypersurfaces are Lie dragged by m^. As consequence of this, if is in E, CmV is 
also in S. 



2.2.3 Eulerian observers 

We can identify the unit timelike vector as the velocity (or the "D-velocity" . . . ) of 
some observer, that we will call the Eulerian observer. The worldlines of these observers are 
orthogonal to the hypersurfaces S^, which means that, for a given t, the hypersurface Tit is 
the set of events that are simultaneous from the point of view of the Eulerian observer. 

We define the acceleration of the Eulerian observer the usual way, 

a^' = n^'Vyit. (2.2.5) 

Let us now list some formula that will be useful for the following sections: 

a^ = V^loga, (2.2.6) 
^fiUo, = -Kai3 - ri/sVa log a, (2.2.7) 
^^m" = -aK"/3 - n^V°Q + n"-^^a. (2.2.8) 



2.2.4 Evolution of 7, 



al3 



From the definition of Lie derivative and equation (2.2.8), one can show 



(2.2.9) 



and 



(2.2.10) 



which means that, for any tensor field T on Sj, its Lie derivative along m is also a tensor 
field on St. 
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2.3 Gauss, Codazzi and Ricci equations 

We still need a way to relate quantities defined on the hypersurface to those defined on 
the manifold Ai; in particular, we would like to have a relation between the D-dimensional 
Riemann curvature tensor, the {D — l)-dimensional Riemann tensor on the hypersurface and 
the extrinsic curvature. Such relations are common in differential geometry — known as the 
equations of Gauss, Codazzi and Ricci — which we now state without proof (see, e.g., [119|). 

2.3.1 Gauss equation 



The starting point is equation (2.1.15). We just need to use equation (2.1.14) to relate VaV' 



with ^ ■ After some algebra we arrive at 



ral''pl'pl''& ''R''a^.u = R^Sap + K\Ksp - pK^s, (2.3.1) 
which is called the Gauss equation. Contracting this equation on 7 and a we get 

T'^al"^ + priori'' ^Rf^.p^ = R^p + KK^^ - K^^^K''^, (2.3.2) 

where we defined K = ^ = K^i (where the equality comes from the fact that K^i, lives on 
E). Taking the trace of this expression we have (noting that K^yK^^ = KijK^^) 



^R + 2'^Rpyn^'rf = R + K^- KijK'K 



(2.3.3) 



2.3.2 Codazzi equation 

We now start with the following equation 



(2.3.4) 



and we project it onto S using (2.1.12). Using equation (2.1.18) and after some algebra we 
arrive at 

l^ral^'prf ''RP^^, = VpK\ - V^K\, (2.3.5) 



which is called the Codazzi equation. Contracting this equation on /? and 7 we have 

^^^c.n" = VaK - V^K^'a. (2.3.6) 



2.3.3 Ricci equation 

We still need one more projection of the Riemann tensor (in fact, the last non-trivial projec- 



tion). Again, we start with equation (2.3.4), but this time we project it only twice onto St 
and one time along n^: 



(2.3.7) 
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Using equations (2.2.7), (2.2.8) and some properties of the Lie derivative we arrive at 



a a 



(2.3.^ 



which is cahed the Ricci equation. We can combine equation (2.3.8) with equation (2.3.2) to 
get 

l^al^p^R^u = —CmK^p - ^Vadpa + R^^p + KK^p - 2Ka^K''p. (2.3.9) 



2.4 Einstein equations 



Our goal now is to write the Einstein equations in an exphcit dynamical system form. Let 
us start by writing the equations themselves in their "traditional" form, 

1 



The alternative form 

where T = g^'^Tf^^, will also be useful to us. 



^TT T, 



T 



D-2 



(2.4.1) 
(2.4.2) 



2.4.1 Decomposition of the stress-energy tensor 



We define 



E = T.^nf^n", 



(2.4.3) 
(2.4.4) 
(2.4.5) 



Sap = T^ul^al'^p, 

which correspond to the matter energy density, the matter momentum density and the matter 
stress density, respectively, as measured by the Eulerian observer. We further define S = 
g'^'S^u = Y^Sij and note that T = S - E. 



2.4.2 Projection of the Einstein equations 



2.4.2.1 Projection onto 



Using equation ( 2.3.9| ) we project equation (2.4.2) onto S^. We get 



Rap + KKap - 2Ka^,K^p + -^^{S - E)^ap - SnSap 



(2.4.6) 

Note that every single term in this equation lives in Sj. Thus, we can restrict the indices to 
spacial ones. 



CmKij = -VjV,a + a 



Rij + KK,j - 2K,kK''j + -^^{S - E)-fij - SnSij 



(2.4.7) 
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2.4.2.2 Projection along n'^ 



This step is easy, we just need to contract equation (2.4.1 ) with n^n*^ and use (2.3.3), yielding 



R + - KijK'^ = IGttE. 
This equation is called the Hamiltonian constraint. 



{2A.i 



2.4.2.3 Mixed projection 



Finally, we need to project equation (2.4.1) once onto Ht and once along n^. Using equa- 

(2.4.9) 



tion (2.3.6 ) we get 



This equation is called the momentum constraint. 



2.5 Choice of coordinates 



Equations (2.4.7) ((Z:* - l)Z:'/2 equations), (2.4.8) (1 equation) and (2.4.9) {D - 1 equations 



contain the same information as equation (2.4.1) (it can be checked that the number of 



independent components is the same: {D — l)D/2 + 1 + {D — 1) = D{D + l)/2). Before 
we can cast these equations in a dynamical system form, however, we have to introduce a 
specific coordinate system, something we have not yet done. 

We will introduce coordinates adapted to the foliation St in the following way |119j . On 

that is varying 
x^~^ is a well 



1 ^2 



each hypersurface we have a coordinate system 
smoothly between neighbouring hypersurfaces, so that rr" 
behaved coordinate system on A^. In this coordinate system 



1 , . . . , 



X 



D~l 



1 9 



We define the shift vector /3 as 



Hfj, = (-a,0, . . . ,0). 



/3 = dt-m, 



(2.5.1) 



(2.5.2) 



or in components, /S'^ = — Note that = 0, so /3 lives on = 0). Using (2.5.2) 

we can also write 



1 



a 



(2.5.3) 



Notice also that 



dfdt = -a' + /3^/3;. = -a^ + P^'jik- 
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We are now able to write the metric components g^iy relative to this coordinate system, 
500 = g^.u (dtT (dtT = dt ■ dt = -a" + 15^ Pu, 

go^ = g^.. {dtf {diT = {m + p)-di = p-di = p^bh = 

9ij = g/iu {diT (djT = Iki {dif {djf = -iij. 
The line element is thus 

g^Ax^dx" = -a^dt^ + (dx' + /3Mt) (dx-'' + P^dt) , (2.5.4) 

or, in matrix form. 



The inverse metric takes the form 



g 



The determinants of g^^ and 7jj are related by 
where 

g = detgf,u, j = det -jij. 



(2.5.5) 



(2.5.6) 



(2.5.7) 



2.6 The PDE system 



Using the properties of the Lie derivative and the definition of shift vector, equation (2.5.2) 
we can write 



^mKij = dtKij — CpKij. 



Equation (2.2.9) can also be put in the form 

{dt - Cp) -fij = -2aKij. 

We now have our full system, which we rewrite here 

{dt - Ci3)jij = -2aKij, 
{dt - Cp) Kij = -ViVja + a 



(2.6.1) 



(2.6.2) 



(2.6.3a) 



■in 



Rij + KK,j - 2K,kK^j + -^^{S - E)jij 



(2.6.3b) 
(2.6.3c) 
(2.6.3d) 
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Note that we can write the covariant derivatives and the Lie derivatives >C/3 in terms of 
partial derivatives of the coordinates x* , and the Ricci tensor Rij and Ricci scalar R in terms 
of and its derivatives in the usual way. This way, assuming that the source terms E, ji, 
Sij are given, we have a second-order non-linear system of PDEs with the unknowns ^ij, Kij, 
a, I3\ 



The above equations (2.6.3) are known in the numerical relativity community as the ADM 



equations, after the work of Arnowitt, Deser and Misner [124j on their Hamiltonian formu- 
lation of general relativity. In this form, however, the equations were in fact first written 
by York |125j (in four spacetime dimensions), and are thus sometimes referred to as the 
ADM- York equations. 

By now we have cast the Einstein equations on an explicit (D — l)-dimensional dynamical 



system form. Note, however, that whereas equations (2.6.3a) and (2.6.3b) are evolution 



equations, equations (2.6.3c) and (2.6.3d) are not. These last two equations constitute 



constraints that the system must satisfy at all times. In particular, these constraints must be 
satisfied at t = 0. So we would now need to specify the relevant initial conditions, satisfying 



equations (2.6.3c) and (2.6.3d), and then evolve them using equations (2.6.3a) and (2.6.3b) 



while making sure that equations (2.6.3c) and (2.6.3d) always hold. 



The question arises: given a specific physical problem (say, a head-on collision of two black 
holes), how does one specify the initial conditions that correspond to the problem we have 
in mind? This is the initial data problem which will be the focus of the next chapter. 



Chapter 3 



Initial data 



On any dynamical system, to perform an evolution one needs to supply the initial condi- 
tions. In our case, this amounts to providing a snapshot of the gravitational fields on some 
hypersurface — the initial data. Then, one evolves this data to neighbouring hypersurfaces 
and so on. 

In general relativity initial data cannot be freely specified. As we have seen in chapter [2| 
not all of Einstein's equations are evolutions equations. We also have a set of constraint 



equations that must be satisfied at all times, the Hamiltonian (2.6.3c) and momentum (2.6.3d) 
constraints. In particular, these equations need to be solved at t = for the {'yij,Kij) that 
represent the physical system we are interested in evolving. We then feed these values to the 
evolution equations themselves. 

In general, this step is far from trivial. There is no unique recipe for the writing of initial 
data corresponding to an arbitrary gravitational system. For some systems, however, — such 
as vacuum spacetimes with moving black holes — recipes do exist. Actually, for the four- 
dimensional case, several methods for constructing initial data for different systems have been 
explored over the years (see |126j for a review). For higher-dimensional systems, however, 
only recently the "standard" way of constructing initial data for moving black holes in the 
vacuum was generalised [127^ I128j . 

In this chapter we will give an overview of the procedure of conformal decomposition first 
introduced by York and Lichnerowicz [129^ \130\ \131\ I132j which rearranges the degrees 
of freedom contained in the three-metric 'jij and extrinsic curvature Kij via a conformal 
transformation and a split of the curvature into trace and traceless part followed by a 
transverse-traceless decomposition of the conformally rescaled traceless extrinsic curvature. 

We will focus specifically on initial data for vacuum spacetimes, generalising the well-known 
Brill-Lindquist [l33] and Bowen-York [IMl [l35] initial data along the lines of [HTl fT28] . 
For alternative procedures to tackle the initial data problem we refer the reader to Cook's 
review [126J, Alcubierre's book |120j . the recent book by Baumgarte &: Shapiro |121j and 
references therein. 
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As in the previous chapter, Greek indices are spacetime indices running from to D — 1; 
Latin indices are spatial indices, running from 1 to D — 1. 



3.1 Conformal decomposition 
3.1.1 Conformal transformations 

We start by recalhng a known general result: given an A^-dimensional manifold with metric 
gfj^i,, if one performs the conformal transformation 

gf,u = (l)ix")g^,u, (3.1.1) 
the Ricci scalars relative to the metrics g^^ and g^^ are related by 

«=i + i-5^Xv.a„0-*^(l-JV)(6-]V). (3.1.2) 

where V is the covariant derivative associated with the conformal metric g^u- 

Let us now consider our case, where we have a (Z?— l)-dimensional spacelike slice with induced 
metric jij and "conformal metric" 'jij. We have N = D — 1, and we make 

We have 



R = 'ip-PR + {2- D)p7p-P-^V''dk^. (3.1.3) 
We further decompose the extrinsic curvature in trace and trace- free parts, 

Kij = Aij + ^^^ij, (3.1.4) 
where K = j'^^ Kij and, by definition, Aij = 0. Defining A^^ = Y'^j^^Aki, we can also write 

K'^ = A'^ + j^^f^ . (3.1.5) 

3.1.1.1 Conformal decomposition of the Hamiltonian and momentum constraint 

Under such a transformation, it is a matter of simple substitution to show that the Hamilto- 



nian constraint equation (2.6.3c) takes the form 



At/; + — — — -R - tA'^Ah - -K^ = 16ttE—^ -, 3.1.6) 

^ p{2-D) p{2-D) ^ p{D-l) p{2-D)' ^ ' 



where A = V'^V 
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With a straightforward calculation we can easily show that 
with q = 2-^^. Thus, we define 

A'^ = iP'^A'^ = iJj^^A'^, 
and we will lower its indices with ^ij, 

We thus have 

ViK'^ = ^p-^ViA'^ + V^K. 

r D-1 



Equation (2.4.9) is then written in the form 



- • D — 2 oD-l ^ . o-D+l 



(3.1.7) 



(3.1.^ 



(3.1.9) 



(3.1.10) 



All we need now is to write equation (3.1.6) as function of A^j, which is very easy. Our system 

D-3 ,D+i 



IS now 



D -2, - D -3 3D-5 
^V' - 777^ 777^ ^ij'^A'^Aij 



A{D-2)^ A{D-2) 



D - 3 D+i 

4(L>-1)^ 



where 



D-2 nP-l 



D-l 



(3.1.11) 
(3.1.12) 



3.2 Initial data for vacuum spacetimes 



Let us now consider the equations (3.1.11) and (3.1.12) for the particular case of vacuum 



solutions (j* = ^ = E). We further impose that the "conformal metric" 7ij is flat (and, thus, 
i? = ) and the "maximum slicing condition", K = {) (to be discussed in section 4.2.1). The 



equations (3.1.11) and (3.1.10) greatly simplify, and we are left with 



diA'^ = 0, 

D - 3 3D-5 . 

^ 4(L» - 2) ^ ^ 



(3.2.1) 
(3.2.2) 



where A = did^ is now the flat space Laplace operator. 
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Here we make a pause to recall that the Schwarzschild-Tangherlini metric in D dimensions is 



r.D-3 



dr + 



dr 



1 L 

rD-3 



+ r^dnl_2, 



(3.2.3) 



where // = '(jyz^xE~2^ ^ being the mass of the black hole and ^at-i = y'(n/2) Siiea, of 
the hypersphere. By performing the coordinate transformation 



r = R(l + 



2 

D-3 



we can write it in isotropic coordinates as 

2 -^^^ V 



ds' 



1 



dt^+(l+ ^ 



We will shortly make use of this geometry. 



[dR^ + R^dn%_2) . (3.2.4) 



3.2.1 Brill-Lindquist initial data 



We now assume that the extrinsic curvature vanishes identically, Kij = 0, a condition that 
holds for time- symmetric initial data. It can be shown |119j that if Kij = and we choose 
coordinates such that a = 1, we have 



^m9o 



0, 



which means that, locally, m'^ is a Killing vector, is also orthogonal to the hypersurface 
T,t=o, and as such this configuration is static. This property only holds locally (on Si=o) and 
we therefore call this configuration momentarily static. 



For Kij = equation (3.2.1) is automatically satisfied, and (3.2.2) reduces to the standard 
D — 1-dimensional flat space Laplace equation, 



AV' = 0. 



We impose the following conditions on ijj 



lim 1^ = 1, 

r— ^oo 



(3.2.5) 
(3.2.6) 



4t 

which is the asymptotic flatness condition (remember that 7jj = ip^'^lij)- 



Let r(j) = \r — X(j)|, where the are arbitrary points in our spacetime. A solution to 
equation (3.2.5) is given by 



N 



1 ^(i) 



(3.2.7) 



where the C(j) are arbitrary constants. Note that equation (3.2.7) obeys the condition (3.2.6) 
The spatial metric takes the form (recall that the conformal metric ■jij is flat) 



N 



-fijdx'dx^ =[l + Yl 



4 

D-3 



i=l ""(i) , 



(dr^ + r-'d^l,^) . 



(3.2.8) 
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This solution is asymptotically flat (by construction), and if we compare this expression 
with (3.2.4), we can identify fi = 4^^^ ^(«)' which is the mass parameter measured in the 
"principal sheet" (anticipating the interpretation). 

We now have to analyse what happens as r — ?■ for a given i. When r — t- r^j) — t- and 
r(j) — >• = \x(^i) — Setting the origin of our coordinate system at r = r(j) (a;(j) = 0) 

we have 



4 

D-3 



1 + 



C(i) 



N 



1 + E 



4 

13-3 



2 ' 



and when r 



(0 



ds^ 



(i) 



4 

D-3 



(0 



1+^ 



(0 



where we defined A( j) = 1 + X] d-3 • With the coordinate transformation r| 



have 



2 

D-3 

(') 



we 



ds^ 



>0 



(0 

{*) D-3 



d^(i)' + ^«'df^D-2;- 



(3.2.9) 



This shows that in this limit the space is also asymptotically flat. Thus, our solution (3.2.8) 
describes a space with -|- 1 asymptotically flat regions. Note that all "lower sheets" are 
separate, i.e., there is no way to travel from one sheet to the other except through the "upper 
sheet" (or "principal sheet"). Equation (3.2.9) shows that each sheet, asymptotically, has a 
Schwarzschild-Tangherlini form, with the mass measured in the ith sheet being given by 



= 4^(i)C(i) = 4 C(i) + ^ 



^D-3 



(3.2.10) 



The observer located on the principal sheet (the {N + l)th sheet) is the only one that sees 



a system of N black holes, with total mass ^adm = fJ-N+i = '^Y2iLi ^(i)^ ^^d already 

mentioned. Thus we identify //(j) = 4C(j) and rewrite our expressions in terms of 



N 



4r 



D-3 ■ 



i=l ^' (i) 

N 



(*)(i) , 



N 



/^ADM 



i=l 



i), 



where r(j) = \r — and = |x(j) — X(^j)\. The points are called punctures. 
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Note that fx / Yli ■ This difference can be attributed to the interaction energy between 
the black holes. It is important to note that //j, as we have defined it, is just a convenient 
label for the mass of the ith black hole (but is not the mass) . The mass of the ith black hole 
as measured on the ith sheet (its "bare mass"), is given by /^(j)!*] 



3.2.2 Bowen-York initial data 



Brill-Lindquist initial data is very useful because it provides us with an analytical solution 
for the constraint equations. However, it also has little physical relevance. Generally, one 
is interested in solutions with black holes that are spinning and moving and as such Brill- 
Lindquist data is clearly not enough. 

In order to have a more general configuration, i.e. one that is not momentarily static, we 
cannot impose Kij = 0. Let us recall that our assumptions are: K = — the maximum slicing 
condition; is flat — the conformal flatness condition; and lim^-^-oo V' = 1 — the asymptotic 
flatness condition. 

We now start by writing A^^ in the form 



A^i = {LXf + i^' 



TT' 



(3.2.11) 



where 



D - 1 



(3.2.12) 



By construction, {LXy^^ij = 0, and we impose ^ijA^r^ = = V jA^r^. We will also restrict 
ourselves to the case ^^rp = 0. The equations (3.2.1) and (3.2.2) take the form 



kx^ + - — ^-d^d^X^ = 0, 
D-l 

D - 3 3D-5 . 

= {Lxy^. 



(3.2.13) 
(3.2.14) 
(3.2.15) 



Thus, we have to solve (3.2.13), plug X^ in (3.2.15) and then solve (3.2.14). We will see that 



even though we will be able to solve (3.2.13) analytically, we generally have to use numerical 



methods to solve (3.2.14). 



To solve (3.2.13) we make the following decomposition ^128j . which introduces functions A 
and Vj, 



X, 



3D - 5, 



D 



(3.2.16) 



* There seems to be some mismatch in the literature as to the definition of "bare mass". Brill and 
Lindquist |133| clearly define it as in our notation, and they even point out that the sum of the bare 
masses is different from the total mass. However, Brandt and Briigmann |135) seem to define bare mass as 
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Equation (3.2.13) then gets the form 



Ay, - x%AV, - 2§-^a, AA - (x^AI4) = 0, 



which is solved if 



AVj = 
AA = 



(3.2.17) 



We have reduced our problem to solving two flat space Laplace equations, which have known 
analytical solutions. In the following we analyse some possible solutions [128J. 



3.2.2.1 Moving black holes 



We start by choosing a solution for the system (3.2.17) of the form 

2tt 



P. 



A = 0. 



(3.2.18) 



iD-2)AD-2 r^-^' 

An stands for the area of the A^-dimensional hypersphere. Pj are constants that, as we shall 
see, will be the linear momentum of the black hole in the j direction. 



For such an ansatz we have, from equation (3.2.16), 

27r 



1 



3L» - 5 



{D -2)AD-2r'^-'^ \D-3 



Pj + {D- 3)n''Pknj 



where nj = 



A'^ 



, and from equations (3.2.15) and (3.2.12) 

Att{D - 1) 1 



n 



pi + n^P' - UkP^^j + {D- 3)n'nip'' 



{D-2)AD-2r^-^ 
The ADM linear momentum is given by |134l 11281 1119j 

^ADM ^ / _ ^d^-^y, 



(3.2.19) 



(3.2.20) 



■in 



(3.2.21) 



where we perform the integration on a hypersphere at infinity; ^/qd y denotes the induced 

dQD-2 (we can 



metric on the hypersphere — using spherical coordinates ^/qd y 
write the induced metric on a hypersphere S of radius r as ds| = qAB^y^dy^ = r^dil|)_2, 
and thus q = det qAB = ('"^)^^^dr2|)_2); is its unit normal vector. 

Reminding ourselves that 

) -D+l 



A'^ = ^^T^^A'^ 

ij ~ ^ 



ij = i + o 



(we are considering K = 0), 
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we see that we can calculate the ADM linear momentum without knowing ip. Plugging ( 3.2.20 ) 



into (13.2.211) we have that i^ADM = p.^ as expected. 



Finally, we note that, as the equation (3.2.1) is linear, we can superimpose solutions of 



the type (3.2.21) corresponding to N Schwarzschild black holes, 

N 



ah 

P{i)^ 



(3.2.22) 



i=l 



where 



{D - 2)Ad-2 



.D-2 
(i) 



(3.2.23) 



where n^^^ = — r^~f^ ^^'-^ parameters P^^^^ correspond to the ADM momentum of the ith. 
black hole when the separation of the holes is very large. 



3.2.2.2 Spinning black holes 



Let us now try a solution of (3.2.17) of the form 



{D - 3)tt Jjk 
{D - 2)Ad-2 ' 



A = 0, 



(3.2.24) 



where Jj^ = —Jkj will be the angular momentum tensor of the black hole. We have 



An-'?. ^ 



and 



A'^ 



^Tl{D - 1) 1 

Ad-2 rD-^ 



(3.2.25) 



(3.2.26) 



The ADM angular momentum is given by (when K = 0) \12S\ I119| 



tADM 



We can check that J^°^ = J, 



""^ Jr—^oo 



Hk- 



For D = 4 we can define the angular momentum vector in the usual way, 

J' = le^'^Jki. 



(3.2.27) 



(3.2.28) 



As in the previous section, we can now superimpose N solutions of the type (3.2.26), 

N 

Af = Y^A%^, (3.2.29) 



1=1 
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where 

^J%) = -^"^^^^ + JS(^i^)h4)) ■ (3.2.30) 

The parameters J^*^^ correspond to the ADM angular momentum of the ith black hole when 
the separation of the holes is very large. 



3.2.2.3 General case 

We can now combine the results from the two previous sections to build a solution of N black 
holes with arbitrary linear momentum and spin, 

N 

A-' = (#(.) + i};.)) , (3.2.31) 

1=1 

\a.b „„j Aab 



where Af^^^ and Af^^^ are given by equations (|3.2.23[) and (|3.2.30[) 



Note: This solution reduces to the Brill-Lindquist momentarily static solution {Kij = 0) 
when P^"^ = = J'^^y For = 1, J"^ ^ and = 0, however, we do not have a slice 
of a Kerr (or, for the higher-dimensional case, Myers-Perry) spacetime. It has actually 
been shown |136) that there is no foliation of the Kerr spacetime that is axisymmetric, 
conformally flat, and reduces smoothly to the Schwarzschild solution in the non-rotating 
limit j^] This means that our Bowen-York solution with J'*^ ^ does represent a rotating 
black hole, but not a stationary one. For the four-dimensional case, when we evolve 
the data, the system emits gravitational radiation and eventually settles down to the 
Kerr solution [137^ 1138] (the higher-dimensional case has not been studied as of yet). 
This spurious gravitational radiation has no desirable physical properties and is often 
referred to as "junk radiation". 



3.2.2.4 Conformal factor 



We still need to solve equation (3.2.14) to get the full initial data, and now there is no hope 



of finding an analytical solution. Let us rewrite the equation we need to solve, 

D - 2, 3D-5 . 



(3.2.32) 



with Aij given by (3.2.31) 



Along the lines of |135j and |128j we write 

"0 = V'BL + 



(3.2.33) 



*For the four-dimensional Kerr solution, but there is no reason to believe that the higher-dimensional case 
is any different. 
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where 



N 

^bl = 1 + EtS^- (3-2.34) 



=1 (t) 



Equation (3.2.32) then takes the form 



D — 3 'I, - 3D-5 

A^x + 4(^_2) ^^^"^^"^ = • (3.2.35) 

For the four-dimensional case, Brandt and Briigmann [135| were able to show the existence 
and uniqueness of solutions for the above equations. Furthermore, the solution for u 
is found on an Euclidean manifold; we do not need to impose inner boundary conditions 
to avoid singularities. Brandt and Briigmann also show that this solution is the "natural" 
generalisation of the Brill-Lindquist initial data, i.e., each puncture represents the infinity of 
another asymptotically flat region of the spacetime and there is no way to travel from one 
sheet to the other except through the "upper" sheet. The higher-dimensional case has not 
been thoroughly studied, but it is believed that the situation is not radically different [128j. 



3.3 Final remarks 

In this chapter we introduced tools for the construction of initial data for higher-dimensional 
numerical relativity. As we mentioned, even though the four-dimensional case has been 
thoroughly studied, the study of initial data for higher-dimensional systems started only very 
recently. As of yet, only Brill-Lindquist and Bowen-York initial data have been generalised, 
but with these two approaches one is already able to construct quite interesting systems for 
vacuum spacetimes. In particular, the Bowen-York approach allows us to write initial data 
for spacetimes with an arbitrary number of moving and spinning black holes. 

For the four-dimensional case there are also powerful computer codes to solve the elliptic 



equation (3.2.35), such as the spectral method presented by Ansorg et al. |139j . 



In upcoming chapters we will present a generalisation of the spectral solver in |139j that 



solves (3.2.35) for black hole binaries in Z) > 5 dimensions with non- vanishing initial boost. 



and preserves the spectral convergence properties observed in four dimensions. 



Chapter 4 

Numerical implementation 



In chapter [2| we have written Einstein's field equations exphcitly in a form (usually referred 



to as ADM equations (2.6.3)) which one could easily give to a computer to evolve. As can be 
seen from this system of equations, though, we are still not quite ready to perform numerical 
evolutions: we still need to say what happens with the variables a (lapse) and /3* (shift). The 
Einstein equations have not imposed any evolution equation for these variables. This reflects 
our coordinate freedom: fixing the lapse function and shift vector is a gauge choice, which 
one could in principle do arbitrarily. In turns out, though, that a good choice is crucial to 
achieve a stable numerical integration. We will in this chapter briefly discuss why this is the 
case and write down the equations we will be using throughout this work. 

It also turns out, as researchers eventually found out empirically in the 1990s when full 
three-dimensional evolutions were attempted using the ADM equations, that this system of 
evolution equations is not well suited to obtain long-term stable numerical simulations. This is 
now known to be due to the fact that the ADM equations are only weakly hyperbolicj^] People 
started experimenting with reformulations of the ADM equations and in 1998 Baumgarte and 
Shapiro — revisiting an earlier formulation based on conformal transformations by Nakamura, 
Oohara and Kojima |140j and Shibata and Nakamura ^141| — showed that this formulation 
behaved much better than ADM for all cases considered |142j . This formulation became 
known as BSSN (Baumgarte, Shapiro, Shibata and Nakamura) and is today the most popular 
scheme used to evolve Einstein's equations. 

It was later realised that indeed the BSSN scheme can be shown to be strongly hyperbolic, as 
opposed to only weakly hyperbolic like in the ADM case, and thus well-posed, e.g. |143l I144j . 

We should also mention that other successful evolution schemes do exit. Most notably, the 
generalised harmonic coordinates approach, e.g. |145j . was successfully used by Pretorius 
in the first ever evolutions of binary black holes through several orbits [18]. Giving a full 
overview of such topics falls outside of the scope of this work. We will in this chapter merely 
introduce the BSSN evolution equations and we refer the interested reader to, e.g., |12mil2l] 



*The ADM equations do allow stable evolutions in spherical symmetry, though (see e.g. |120) V 
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for comprehensive overviews. 

We close this chapter by introducing the numerical code itself used for all the simulations to 
be presented. 

In this chapter, we restrict ourselves to the four-dimensional case (for completeness, we 
present in appendix 4. A the higher-dimensional BSSN scheme), and therefore spatial (Latin) 
indices are here restricted to i = 1, 2, 3. 



4.1 BSSN formulation 

As we have just mentioned, if we were to try and evolve Einstein's equations in the ADM 



formulation (2.6.3) (supplemented with the gauge conditions we will introduce in the next 
section) we would find out that the system is severely unstable. In this section we recast the 
evolution equations in the BSSN form, which allows for stable numerical evolutions. 

We start by performing a conformal decomposition of the spatial metric 'jij in the following 
way (observe it is going to be a different decomposition from the one performed in the study 
of the initial data) 

= Xlij ■ (4.1.1) 

The conformal factor x can in principle be freely prescribed. In the BSSN scheme, one 
imposes that the determinant of the conformal metric be equal to the determinant of the flat 
metric rjij, 

By construction, we have 

det% = 7]. (4.1.3) 

Since we will stick to Cartesian coordinate systems throughout this work, we will always have 
r] = 1 = det'jij, which makes x a scalar density with weight —2/3. 



-1/3 



(4.1.2) 



Just like in equation (3.1.4), we decompose the extrinsic curvature Kij into trace and trace- 
free parts and apply the conformal transformation we used for the metric to the traceless 
part, 

Kij = + f • (4-1-4) 

Let us now find evolution equations for the variables we have introduced {x:7ij:K,Aij). 



(4.1.1 


) and ( 


4.1.4) into ( 


2.6.3a 


) and ( 


2.6.3b 



and 



dtX = l3''dkX + ^x{aK-dkPn, 



dtK = P^dkK - V^dua + a (^A'^ Atj + \k^ ] + 4vra(^ + S), 



(4.1.5) 



(4.1.6) 
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where in this last equation we have already used the constraint (2.6.3c) to eliminate the Ricci 



scalar. Substituting back we can compute the remaining evolution equations, which take the 
form 

dtjij = P^dk% + 2%^idjjP'' - ^%dkl3'' - 2aAij, (4.1.7) 

and 



+ a {KAij - 2Ai''Akj^ - Svra (xSij - 



(4.1.8) 



where '^^ denotes the trace-free part, e.g., RJ^^ = Rij — ^JijR- 
We further need to decompose the Ricci tensor in two parts. 



Rij - Rij + Rfj 



(4.1.9) 



where Rf^ only depends on x Rij is the Ricci tensor associated with the metric ^ij. 
This term contains mixed second derivatives of the metric, something that is undesirable. 
For stable numerical integration, the following "conformal connection" variable was intro- 
duced dm [Ti2] 

r = y'=f;., = -a,f^-. (4.1.10) 

In terms of this conformal connection, the conformal Ricci tensor then takes the form 



Ri 



(4.1.11) 



As we can see, the first term in this expression, which involves a Laplacian, is the only explicit 
second order derivative operator acting on ^ij. All the mixed second derivatives have been 
absorbed in derivatives of F*. Since the BSSN scheme considers P to be an independent 



variable, we need an evolution equation for it. Acting on (4.1.10) and interchanging the time 
and space derivatives we get 



dtV = -dj ( p'^dkf^ - 2f^idkp'^ + -f^dup'' + 2aA'^ 



(4.1.12) 



We further use the momentum constraint (2.6.3d) to do away with the divergence of the 
extrinsic curvature and obtain 



- ^af^djK - A'^ {3ax-^d,x + 2d,a) - IG^x^ V- 



(4.1.13) 
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The full system of evolution equations is then 

dtlij = P'^dk'rij + 27fe(i5,)/3'= - ^^ijdkP'' - 2aAij, (4.1.14a) 

dtX = P''dkX + lxiaK-dkl3''), (4.1.14b) 

dtAij = p'^dkAij + 2ifc(ia,)/3'= - '^Aijdkl3'' + X {aR^J - ViS^a)^^ 

+ a (if i,, - 2Ai^Akj) - Svra (^xSij - , (4.1.14c) 

dtK = li^dkK - V^OfcO + a {a'^ Aij + ^K^^ + 47ra(S + 5), (4.1. 14d) 

- \af^djK - A'^ {Sax'^djX + 2dja) - 16^ax" V , (4.1. 14e) 
where = + -R^- 



R 



X 



\x-^ (drd.x - dkXT\^ - \x-^d,xd,X - \l^■JX-^ dkX^'' (4.1.15) 



+ li^n^'x-' (dkdix - • 

Source terms are determined by 

Sij = 7"a'^jTa[s , s = Y^Sij . 



(4.1.16) 



The above system of evolution equations (4.1.14) is known as the BSSN evolution scheme 



and has proven to be extremely robust for numerical evolutions of Einstein's field equations. 
Numerous other schemes do exist; most, however, offer no substantial advantage over BSSN, 
which has remained extremely popular. We will use BSSN for all our numerical evolutions. 



4.2 Gauge conditions 

We now turn our attention to the gauge conditions. The first question to ask is: what is a 
"good" choice for a and /?*? An obvious first choice, also the simplest possible, is to impose 
the so-called geodesic slicing (also known as Gaussian normal coordinates), 

a = l, I3' = 0. (4.2.1) 

This choice was in fact used in the pioneering work by Hahn and Lindquist [8]; it is now 
known, however, that it is actually a very bad choice for long-term evolutions. We can 
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intuitively understand why this is the case. First we note from equation (2.2.6) that the 
Eulerian observers have zero acceleration and thus follow geodesies (hence the name of this 
slicing). In the presence of black holes (or other gravitational sources), geodesies tend to 
focus. Coordinate observers will then collide with each other, consequently forming coordinate 
singularities and crashing the numerical evolution. We thus need better gauge choices. It 
falls outside the scope of this work to give an overview on the merits and disadvantages of 
the different conditions that have been proposed throughout the years. We will simply state 
and motivate the conditions we will be using. 



4.2.1 1+log slicing 

A famous choice for the lapse function is known as maximal slicing, which corresponds to 
imposing that the trace of the extrinsic curvature vanishes throughout the evolution, 

K = 0. (4.2.2) 

A nice property of this condition is its singularity avoidance. We can see this by taking 



the trace of equation (2.1.18), which with (4.2.2) implies VnU^ = 0, an incompressibility 



condition on the velocity field of the Eulerian observers. This prevents the observers from 
converging and the subsequent appearance of a coordinate singularity, as in the geodesic 
slicing case. Such a property is very much desirable, making maximal slicing an attractive 
choice. There is an enormous disadvantage, however, which is the need to solve an elliptic 



equation at every time step during the numerical evolution in order to ensure (4.2.2). We 
therefore would like to have conditions that mimic this property of maximal slicing, yet with 
a hyperbolic character. 

Such a choice is the so-called 1+log slicing, 

{dt - Cp) a = -2aK , (4.2.3) 

which, being a hyperbolic equation, is trivial to implement numerically, has been shown to be 
extremely robust and mimics the singularity avoidance properties of maximal slicing |146j . 



This condition gets its name from the fact that, when imposing /3* = 0, equation (4.2.3) can 
be integrated to give 

a = 1 + log 7, (4.2.4) 

where we recall that 7 = det jij . 
4.2.2 Gamma driver 

Having chosen a condition for the lapse function, it remains then to say what happens to the 
shift. 

A possible choice, known as the Gamma freezing condition, is the following 

dtr = 0. (4.2.5) 



CHAPTER 4. NUMERICAL IMPLEMENTATION 



50 



Using equation (4.1.14e), we can write the above condition as an elliptic equation for /3*. 
This condition is related to the "minimal distortion" shift condition |147| . which attempts 
to choose coordinates such that the time derivative of the 3-metric dt'jij is minimised. The 
disadvantage is once again the need to solve an elliptic equation at each time step. 

Researchers have therefore proposed alternative conditions, using parabolic or hyperbolic 
equations, that mimic the minimal distortion condition with good approximation. The 
following choice (and variations thereof) is now extremely popular 

{dt-Cp)f3' = r-Vpf3\ (4.2.6) 

where r/^ is a function of spacetime. This is known as the Gamma driver condition [148J. 

Use of these gauge choices proved crucial for the 2005 breakthroughs using the moving 
puncture technique |20| [T9] . 



4.3 Numerical code 



Having chosen a set of evolution equations (4.1.14), gauge conditions (4.2.3), (4.2.6) and 



prescriptions for setting initial data (see chapter [3| , it remains then to assemble everything 
on a numerical code. Such a task is far from trivial. One of the main reasons is that the 
presence of very different scales in the spacetimes that are usually evolved requires the use of 
mesh refinement. The problem is further complicated by the need to use parallel computing 
and to store large amounts of data. 

All numerical simulations that will be presented in subsequent chapters have been performed 
by adapting the Lean code |149j . initially designed for 3 + 1 vacuum spacetimes by U. Sper- 
hake. Lean is based on the Cactus computational toolkit [150], it employs the BSSN 
formulation of the Einstein equations |14H I142| (with fourth order discretisation in the spatial 
derivatives) with the moving puncture method |20tll9j. uses the Carpet package for Berger- 
Oliger mesh refinement [151^ 1152] . the spectral solver described in |139j for 3 + 1 initial data 



and Thornburg's AHFinderDirect |153t 1154] for horizon finding (see section 5.2) 



For a given numerical simulation our numerical grid will consist of two types of cubic 
refinement levels: n outer levels centred on the origin (remaining stationary throughout 
the simulation), and m inner levels centred around each black hole (and following these as 
the simulation progresses). The following notation (which we will make frequent use of) 

{(256, 128, 74, 24, 12, 6) x (1.5, 0.75), h = 1/48} 

specifies a grid with six fixed outer components of "radius" 256, 128, 74, 24, 12 and 6 
respectively and two refinement levels with two components each with radius 1.5 and 0.75 
centred around either black hole. The resolution is h = 1/48 on the finest level and 
successively decreases to 2'^/48 = 8/3 on the outermost level. Further details about Lean 
may be found in |149j . 
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4. A D-dimensional BSSN equations 

For completeness, we here write the full L'-dimensional BSSN equations, as first written 
in [27] . These equations can be derived in a procedure entirely analogous to the one outlined 
in section |4?T1 



daij = P'^dkiij + 27fe(^c>,)/3'= - -^—^l^j^kP'' - 2aAij, (4.A.la) 

dtX = li''dkX+^-^x{aK-dkP''), (4.A.lb) 
dtAij = p^dkAij + 2ifc(,a,)/3'= - -^—^AijdkP'' + X {aR^j - Vidjaf' 

+ a(^K Aij - 2Ai^Akj^ - Sira (^xSij - jyzi'^i^j , (4.A.lc) 



dtK = (5^dkK - V'^dta + a ( A'^Aij + - — -K^ ] + ^^^a [{D - 3)E + 5] , (4.A.ld) 



D-1 J D-2 

- 2^^af - A^^ (^{D - 1)"^ + 29ia) - 167rax"V, (4.A.le) 

where denotes the trace-free part and the Ricci tensor Rij is further split into Ri 
Rij + where 



- - . , , ■\itV~i rwa 



^ zM Pi^ , z, fifc pi_ z, a zhl I ^ fifc z. f^k 



Rfj = ^^X-' (^^^,x - dkxrf,) - ^^X-'d,xd,x - hjX-'dkXr" (4.A.2) 



+ \l^3l'''x-' (dkdix - ^-^^-Y^X-'dkxdiX 



Equations (4.1.14) can be recovered with D = A. 



Chapter 5 



Wave extraction and horizon 
finding 

We have thus far covered, essentially, all the main tools necessary to successfully evolve 
Einstein's equations on a computer. Assuming then that we can specify some arbitrary initial 
configuration and evolve it for as long as we like, we are still faced with the most important 
task: how to extract the relevant physical information. Recalling that the coordinate system 
used throughout the evolution is designed to be well suited to the numerical evolution and 
not for human-readability, we easily convince ourselves that it is not trivial to read physical 
information from the numerical output. For this purpose, tools were developed to enable 
the extraction of the gravitational wave information from a numerical simulation and, when 
dealing with black hole spacetimes, information about the black hole's horizon. 

In this chapter we will briefly outline the two main methods of extracting gravitational wave 
information and the corresponding waveforms: the gauge invariant formalism of Kodama and 
Ishibashi |155l I156| — itself a generalisation to higher-dimensional spacetimes of the Regge- 
Wheeler-Zerilli formalism |157| 1158] , later put in a gauge-invariant form by Moncrief |159] — 
and the Newman-Penrose formalism |160) . We will also mention the very basics regarding 
finding (apparent) black hole horizons. 



5.1 Wave extraction 

Gravitational waves are ripples in the shape of spacetime that propagate information at finite 
speed, just as water waves are ripples in the shape of an ocean's surface. They are one of 
the most important predictions of general relativity. These waves have never been directly 
detected; there is, however, strong indirect evidence for its existence since the discovery of 
the famous binary pulsar PSR 1913+16 (also known as the Hulse- Taylor binary pulsar after 
its discoverers |161| ). whose orbital period change is consistent with the general relativistic 
prediction for energy loss via gravitational wave emission. Other systems have since been 
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uncovered, allowing for even more stringent tests, e.g. jl62t I163| . 

These waves are generated by dynamical gravitational fields — roughly speaking, accelerated 
bodies in non-spherically symmetric motion will emit gravitational wave^ As they propagate 
throughout spacetime, they carry with them information about the physical properties of the 
system that produced them. By measuring them with gravitational wave detectors — such as 
the already mentioned LIGO, Virgo, GEO and TAMA — we can have a brand new window 
opening up to the universe. The likelihood of such detections is greatly enhanced if one 
can use theoretical gravitational wave signals coming from possible astrophysical sources 
as templates. Our task here is to briefly introduce the techniques used to generate such 
templates from numerical simulations. 

Before we begin, let us make one last comment. We have mentioned energy carried away by 
gravitational radiation, but as we know there is no notion of local energy of a gravitational 
field, so some care has to be taken here. The usual procedure is to write a stress-energy 
tensor for the metric fluctuations that is second-order in said fluctuations (the same way 
that the stress-energy tensor associated with a scalar or electromagnetic field is second order 
in the fields). Modulo some subtleties, such a quantity can be constructed, and meaningful 
quantities can be extracted from it. We will not be giving details on its derivation or the 
subtleties involved (see, e.g., |1641 1165| I166j ). but merely present the relevant formulae that 
will be of use to us. 

We will start by recalling known four-dimensional results that will be of use. In the weak-field 
limit, we can write the metric tensor as the Minkowski metric plus perturbations, 

g^iu = V^u + h^u , \h^u\<^'^- (5.1.1) 
To first order in /i^i,, the Riemann tensor is 

- {di3dfj,hau + dadyhp^ - dpdyhaf, - dad^h^y) . (5.1.2) 
It is useful to introduce the usual trace reversed perturbation 

hf,u = h^u - ^r]^y . (5.1.3) 

Imposing the Lorenz gauge 

d^h'''' = 0, (5.1.4) 

the linearised field equations reduce to 

□V = -16^r^^, (5.1.5) 

where □ is the d'Alembertian operator in flat space. In vacuum we get the usual wave 
equation 

□V = 0. (5.1.6) 

*For a system to emit gravitational waves the third time derivative of its quadrupole moment has to be 
non-zero. 
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Since the Lorenz gauge does not completely fix our degrees of freedom, we can further impose 
the transverse-traceless (TT) gauge 



0. 



(5.1.7) 



where, to simplify, we can use a Cartesian coordinate system rj^i, = diag(— 1, 1, 1, 1) and u'^ is a 
unit timelike vector. The second equation reflects the fact that there is no propagating scalar 
mode in general relativity. Note that in this gauge h^i, = h^^. With the constraints (5.1.4) 



and (5.1.7) we are left with two degrees of freedom (in four dimensions; generically, in D- 
dimensions, we have D{D — 3)/2 degrees of freedom). We can write the plane- wave solution 



of equation (5.1.6) subject to the constraints (5.1.7) in the usual form 



h 



TT 



(5.1.8) 



taking = (1, 0, 0, 1) and where 



A 



(0 








o\ 





h+ 


hx 








hx 















0/ 



(5.1.9) 



/i+ and hx are the two independent polarisations of the gravitational wave, known as "plus" 
and "cross" polarisations. 

It can be shown (e.g. |166tll23"tll21j ) that the outgoing energy flux carried by the gravitational 
radiation is given by 



dt 



lim 

r-5>oo \Qtt 



h\ + hl 



(5.1.10) 



To derive this formula, we need to expand Einstein's equations up to second order perturba- 
tions. Terms that are quadratic in the first order perturbations of the metric, after suitable 
averaging, can then be viewed as sources, constituting an effective stress-energy tensor for 
gravitational waves. This stress-energy tensor can then be used to compute energy and 
momentum carried away by the gravitational radiation. 



5.1.1 Newman-Penrose formalism 

We now briefly describe the Newman-Penrose formalism. This formalism (also known as 
the spin- coefficient formalism) introduced by Newman and Penrose in 1962 [160] is an 
alternative way to write Einstein's equations which has proven to be extremely useful in many 
situations in general relativity, such as in searches of exact solutions, black hole perturbation 
theory and studies of gravitational radiation. There is a whole literature devoted to this 
formalism. For its application in numerical simulations, we mention for instance the books 
by Alcubierre |120j . Baumgarte &: Shapiro |121j and references therein. Here we will only 
state the basic equations that we will need and refer the reader to relevant publications where 
appropriate. 
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In this section we restrict ourselves to four-dimensional spacetimes since this formalism has 
not been generalised to higher dimensions j^] 

This formalism starts with introducing a null complex tetrad k,m,rri} satisfying 

— I ■ k = 1 = m ■ fh , (5.1.11) 

where rh is the complex conjugate of m and all other inner products vanish. 

We further note that the four-dimensional Riemann tensor '^^^i?"^^^ has 20 independent 
components. Its trace, the Ricci tensor ^^^Rap has 10. The remaining degrees of freedom 
are encoded in the Weyl tensor *-^^Ca/375, defined as 

= - (^k[,(%5]/3 + + \^%a[,^'W^R ■ (5.1.12) 

The Newman- Penrose formalism encodes these degrees of freedom in a set of complex scalars, 
often called Newman-Penrose scalars. The ten independent components of the Weyl tensor 
are encoded in the five complex scalars ^'oj . . . ; ^4 (often also called Weyl scalars)^ All of 
these scalars are formed by contracting the Weyl tensor (and the Ricci) with the complex 
null tetrad. Since there is no unique choice for a null tetrad satisfying (5.1.11), the choice of 
this tetrad will affect the Weyl scalars and their physical interpretation. 

For a class of such tetrads, the so-called quasi- Kinnersley frames, and ^'3 both vanish, 
and we can interpret ^'o and ^4 as measures of the incoming and outgoing gravitational 
radiation, whereas ^'2 can be interpreted as the "Coulombic" part. ^'4 and ^0 defined 
ai 

'^o = ^^'^Ca.p^sl^m^Prn^ ■ (5.1.13) 
4 = ^^^C^p^sk'^rfil^k^'m^ . (5.1.14) 

Since (for the suitable tetrad we mentioned above) the latter quantity encodes the outgoing 
gravitational wave signal, this will be of particular use to us (^0 will also be of use in 
section 



7.2). 



In practice, we construct /, k and m from an orthonormal triad 6^,6^,6? orthogonal to the 



unit timelike vector e£ 



1 

1 



I = —f=. {ei + ef 



k = ^{ei-er), (5.1.15) 



*A related formalism also based on spin-coefScients, the Geroch-Held-Penrose (GHP) [167] formalism, has 
been extended to higher dimensions [168| . 

^The ten independent components of the Ricci tensor are analogously written in terms of four scalars and 
three complex scalars, but we will never make use of these quantities in this work. 

'''Different sign conventions exist in the literature. 
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We refer the reader to |169j for a review of the formahsm; here we merely note that asymp- 
toticaUy the triad vectors e^, e^, behave as the unit radial, polar and azimuthal vectors. 

Having chosen our tetrad, we can now compute an explicit expression for ^'4 using the 
definition (5.1.14). In the TT gauge, this can be shown to be, for outgoing waves \120\ 1121] 

^'0 = 0, (5.1.16) 

^4 = + i/ix , (5.1.17) 

whereas for ingoing waves, we have instead 

^-0 = /i+ - i/ix , (5.1.18) 
^4 = 0, (5.1.19) 

where ' denotes a time derivative and h+ and hx are the amplitudes of the plus and cross 
polarisation of the gravitational wave (5.1.8), (5.1.9). Herein lies the usefulness of the ^'4 
scalar. 

It is useful to perform a multipolar decomposition by projecting ^4 onto spherical harmonics 
of spin weight s = —2 (cf., e.g., appendix D of [12U]): 

^4(t,^,0) = j;V''™(i)>^J(^,<^) ■ (5.1.20) 

Lm 



In terms of these multipoles, the radiated flux is given by the expressions \160\ I170j 



dE'cw r 



dt 



lim 
r^cx) 167r 



E 

Lm 



(5.1.21) 



In Einstein-Maxwell theory, the right-hand-side of Einstein's equations reads 



T, 



1 



Efi ^g^yF F\fj 



(5.1.22) 



where F^^'^ is the Maxwell-Faraday tensor. In such cases, we can analogously extract the 
electromagnetic wave signal in the form of the scalar functions, and $2 [160^ 1170] . deflned 
as 



$2 = F^^m'^k'' . 



For outgoing waves at inflnity, these quantities behave as 



~ ^ {Ef + iBf 



$2 



E2 - lEi 



(5.1.23) 
(5.1.24) 

(5.1.25) 



Again, it is useful to perform a multipolar decomposition by projecting $1 and $2 onto 
spherical harmonics of spin weight and —1 respectively: 



ci>i(t,^,0) = j;</.f(t)y,o„(0,<A) , 

l,m 



(5.1.26) 
(5.1.27) 



Lm 
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In terms of these multipoles, the radiated flux is given by the expressions [16^ I17U| 

.2 , 2 



FEM = ^=lim^j:Ur(t) 

Lm 



dt 



(5.1.28) 



We see from (5.1.25) that ^2 encodes the radiative modes. 



5.1.2 Kodama-Ishibashi 

A different approach to extract gravitational wave perturbations is that of the gauge-invariant 
Moncrief formahsm |159j . This has been generalised to higher dimensions by Kodama and 
Ishibashi (KI) |155l I156j , and we will review this approach in the following. 

In the KI formalism, we start by writing the metric element as a background metric plus a 
perturbation 

gAB = g% + 6gAB. (5.1.29) 
The background spacetime has the form 

ds2(0) = gf^dx^dx^ = gi°^dx"dx^ + r'^dnD-2 = ^i^dx'^dx'' + r^^J^^dx^x^ , (5.1.30) 

where the x^ coordinates refer to the whole spacetime {A = 0, . . . , D — 1), x"" = t,r and il^^ 
is the metric on the unit {D — 2)-sphere S^~^ . 

The procedure now is to expand the metric perturbations qab into harmonic functions. These 
exist in three flavours — scalar, vector and tensor harmonics. Metric perturbations can then 
be written in terms of gauge invariant quantities |155j : 

Tensor harmonics T^^; satisfy 

(A + fc2)T,5 = 0, (5.1.31) 

with the properties 

T\ = 0, T\^ = 0. (5.1.32) 

where A is the Laplace-Beltrami operator on S^^^ and -a denotes the covariant deriva- 
tive with respect to the metric 0.^^ on the sphere. Note that we will omit the index 
labelling the harmonic throughout this discussion. 

For tensor-type perturbations, the metric perturbations Sqab = h^B are expanded in 
the following way 

hab = 0, haa = 0, h-^i = 2r^ HtT -^i , (5.1.33) 

where Ht = Hxit, r) (again, we leave the harmonic labels implicit), and note that there 
is sum over the indices of the harmonics in this expression. 
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Vector harmonics Vq satisfy 

(A + A:2)Va = 0, (5.1.34) 

V^:a = 0. (5.1.35) 

from this, vector-type harmonic tensors can further be defined 

We expand the vector-type perturbations as 

hab = 0, Ka = rfaYa, h-,i = 2r'' HtY -,f, , (5.1.37) 
where fa = fa{t,r). 
Scalar harmonics S satisfy 

''a + A;2)s = 0, (5.1.38) 
from which we can build scalar-type vector harmonics 

§a = -^S:a, (5.1.39) 

and scalar-type harmonic tensors 
We expand scalar-type perturbations as 

hab = fabS , ha-a = rfaSa , h-,l = 2r^ [Hl^-^i^ + Ht\i) , (5.1.41) 

where fab = fab{t,r). 

For Z > 1, the metric perturbations can be expressed in terms of the following gauge- 
invariant variables |155j 

F = Hl + ^^Ht + ^ X.rl'^ , 

D-2 r (5.1.42) 

Fab = fab + + X^ia , 

where we have defined 

= ^ [fa + ^i^TIa) , (5-1.43) 

and we denote the covariant derivative with respect to the metric g'^^^^ with a subscript 

\a- 
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Impressively, a master function <^> can be defined that, from the perturbed Einstein equations, 
can be shown to obey the simple wave equation (T56j 

(□-y(r))$ = 0, (5.1.44) 

where □ is the d'Alembertian operator with respect to g^^J , and the form of the potential 
V{r) depends on whether one is considering scalar, vector or tensor perturbations. 

The master function <I> is specially useful since it encodes the gravitational waveform; the 
energy emitted via gravitational radiation can also be computed quite effortlessly. Writing 
the index I explicitly, the energy flux in each Z-multipole is |171j 



^ = ^^k\k' - D + 2){^^,^' 
dt 327rL>-2 ^ '* 



k'^ik^ - D + 2){<^>\f . (5.1.45) 



The total energy emitted in the process is then 



1=2 



5.2 Horizon finding 



When evolving black hole spacetimes, besides the wave extraction tools, physical information 
can also be read from its horizon properties. A black hole is a region of spacetime from which 
no future directed null geodesic can reach an outside observer. Its surface, the event horizon, 
acts therefore as a one-way membrane. In asymptotically flat spacetimes, the event horizon 
can be defined as the boundary of the causal past of future null infinity. It is thus as global 
concept, requiring information from the whole spacetime to be located. From the point of 
view of a numerical evolution, this is not very useful since one would like to know about the 
location of the black hole as the simulation progresses. 

A more useful concept in this regard is that of the apparent horizon. It is defined as the 
outermost marginally trapped surface on a given spatial hypersurface — a closed surface on 
which the expansion of (outgoing) null geodesies vanishes. The apparent horizon is a local 
concept, depending only on information present on the given hypersurface, making it an ideal 
diagnostic tool for numerical evolutions. 

Given a spatial section S of a spacetime with 3-metric "jij and extrinsic curvature Kij, the 
expansion of null geodesies can be shown to be 

G± = ±Vis' + Kijs's^ - K (5.2.1) 

where V is the covariant derivative with respect to the 3-metric ^ij and s' is the spatial 
normal to the apparent horizon surface within S. 0+ is the expansion of the outgoing null 
geodesies, 0_ the expansion of the ingoing ones. The (black hole) apparent horizon is then 
defined by the following equation 



Vis' + Kus's^ - K = {). 



(5.2.2) 
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General purpose tools exist to solve this equation during numerical evolutions; for an overview 
see e.g. ^172j and references therein. 

Finally, we emphasise that apparent horizons are slicing dependent. It is possible, for instance, 
to foliate the Schwarzschild spacetime in such a way that there is no apparent horizon |,173] 
(the event horizon, being a global quantity, is an intrinsic property of the geometry and is 
thus always present). The presence of an apparent horizon, however, does imply the existence 
of a section of an event horizon exterior to it (assuming cosmic censorship and R^^k^^k'^ > 
for ah null k^' [174]). 



Chapter 6 



Higher-dimensional numerical 
relativity 



As mentioned in the Introduction, the abihty to perform fully non-linear numerical evolutions 
of Einstein's field equations in higher-dimensional scenarios has tremendous potential to 
answer fundamental questions in physics, with possible applications including studies of the 
AdS / CFT duality, explorations of TeV-gravity scenarios and the study of higher-dimensional 
black hole solutions. 

Numerical relativity in higher dimensions has only recently started being explored, with 
pioneering works including those in [92l [971 1175^ [571 1176j . In this chapter, we will describe 
the approach of [175^ 1176] . 

The formalism we will present allows us to consider two classes of models, which are gen- 
eralisations of axial symmetry to higher dimensional spacetimes: a D > 5 dimensional 
vacuum spacetime with an SO{D — 2) isometry group, and a D > Q dimensional vacuum 
spacetime with an SO{D — 3) isometry group. The former class allows studies of head-on 
collisions of non-spinning black holes. The latter class allows to model black hole collisions 
with impact parameter and with spinning black holes, as long as all the dynamics take 
place on a single plane. This class includes the most interesting physical configurations 
relevant to accelerator — and cosmic ray — physics (in the context of TeV-scale gravity), and 
to the theoretical properties of higher-dimensional black objects (such as stability and phase 
diagrams) . 



In section |6.1[ we introduce a general dimensional reduction procedure from D-dimensional 
vacuum general relativity to a lower dimension model; in section |6.2| we specialise the 
equations obtained to the case where the Z)-dimensional spacetime has an S0{D — 2) isometry 
group, perform the 3+1 splitting of space and time and write down a system of evolution 
equations; in section |6.3| we outline the construction of relevant initial data, following the 



approach of |177] : in section 6.4 we discuss some code tests, introduce a wave extraction 



procedure and present results. We end this chapter with a discussion in 6.5 
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We note that in this chapter, due to the necessity of introducing multiple covariant derivatives, 
we shall explain the notation as we go along. 

6.1 Dimensional reduction 

The starting point of our formalism is a dimensional reduction from D-dimensional general 
relativity in vacuum to a lower dimensional model. 

The isometry group of a Schwarzschild (or, for D > 4, Tangherlini |178j ) black hole is 
SO{D — 1) X M. For a head-on collision of two non-rotating black holes, the isometry is 
further reduced to SO{D — 2): indeed, neither the time direction nor the direction of the 
collision correspond to symmetries, but a rotation of the remaining D — 2 spatial directions 
leaves the spacetime invariant. 

One can take advantage of this symmetry to reduce the spacetime dimensionality. This can 
be accomplished by writing Einstein's equations in D dimensions in a coordinate system 
which makes the symmetry manifest, allowing for a lower dimensional interpretation of the 
D-dimensional Einstein's equations (in the spirit of the Kaluza-Klein reduction). We remark, 
however, that we do not perform a compactification; rather, we perform a dimensional 
reduction by isometry, as first proposed by Geroch |179j . The extra dimensions manifest 
themselves in the lower dimensionality as a source of Einstein's equations, defined on the 
lower dimensional manifold. 

In the original proposal of Geroch ||179j the symmetry space was SO{2). This approach 
has been applied to numerical relativity, see for instance [180^ I18H I182j ; a five dimensional 
extension, with the same symmetry space, has been derived in |183j . A generalisation to 
coset manifolds (like the sphere 5") was given by Cho in |184[ 1185] , but in these papers the 
complete form of Einstein's equations was not presented. 

Following the approach by Cho |184j . we will start by deriving the general equations obtained 
doing a dimensional reduction by isometry. We will afterwards focus on the isometry group 
of the S*" sphere and present the equations obtained with a dimensional reduction to four 
dimensions, as well as their numerical implementation. 

6.1.1 General formalism 

The most general D-dimensional metric gAB, ^ = 0, . . . ,d — 1, . . . ,{D — 1), can be written 
in the following form (in the coordinate basis = {d^, di)) 



ds 



QAsdx'^dx 



.B 



i- 




where n = 0, . . . ,d — l and i = d, . . . , D — 1. Kisa scale parameter and e a coupling constant. 
This metric is fully general and not an ansatz. 
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Assume that g^B admits an m-dimensional isometry G, generated by m Killing vector fields 
which we express as (assuming D — d is large enough) 

= KQl. (6.1.2) 

a = 1, . . . , m = dim G. The Killing vector fields form the Lie algebra of G, satisfying 



Defining 



^^a9AB = 0, 

1 

[Caj'^fe] = ~/ ab^c- 



[iaAi\^Fi.d-j = -[(kKi)d-^, 



and the "dual" form to the Killing fields by 



5.1.3) 
5.1.4) 



we can derive, from (6.1.3), 



where F t 



(hg-fk - F ijOik + F ikg-ji, 
digt,u = 0, 



(6.1.5) 



Our goal is to compute the Ricci tensor of metric (6.1.1), which is more easily done in a 
non-coordinate basis. Details of the computation can be found in appendix 6. A here we 
mention only the final result. 

We first define the "covariant derivatives" and as 



V^T'- 
ViT 



kfj, 



T~) rpioi _ I -772 j~nla_ 'T'l _'~pioi_ i "pQ '~pi^_ "pA rpia_ 
^(y^ kii^ al^ kii ak^ Ifi^^ Ao-J k^l ^ Mo--i fc^' 



j-^ k^i = '-'j'^ fc/x J- Ij-^ kfi~ kj-'- Ifj, 



(6.1.6) 

(6.1.7) 



where 



-eKF\-Bl, 



J""- = 



-bkG^ 



-eK{d^B',-d,B'^ + eKt'--B;Bi. , 



ft__ — -pi 7?*-- 

jk ~ jk ^ kjt 



and both connections are metric, 

= dj^gij — T^jkQYj — ^''jkdii = 0. 



(6.1.^ 
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Note however that 



The Ricci tensor of (6.1.1) is (see appendix 6. A) 



.1.9) 



R^ = bkR^bI + ^g^^V^ {g-fkJ'K,) + \g'''VpgMg^^r^x,mm + ^^k {i'^V^gu 



(/V^ff^^) = R-^, (6.1.10) 



R^u = R^u + 2eKBlR,-)j - e'K'R^-Bl^Bi - -g'^^gij^ x^T^ au - -V. ( g^^V^^IJ 



Ig^'i^'yf.g-fk'^ug-fl- )^-k^~\.. (e.i.n) 



and 



R = R + R- \g-,jg-'^g^''F\,FKu - [g^^V ,g-,^ - \g^' g^'V ^ gifJ ,913 

- -/'g^^V^^g-kf^^.gi3- (6-1-12) 
These are the expressions we were looking for. Equivalent forms can be found in |183| I184j . 



6.1.2 Examples 
6.1.2.1 51 

As a first (trivial) exercise, we can reproduce the standard Kaluza-KIein expressions. Re- 
member that the Kaluza-Klein metric has the form 

ds2 = g^^dx^dx'^ + e^'t' (drr^ + A^dx^'f . (6.1.13) 

We can easily recover this case from our formalism by making d = A, D = 5, gij ^ e'^'^, 
g^^ e~'^^, enB^^ A^, and T\u -F^u = - {d^A^ - duA^) (cf. equations ( |6.1.8[ )). 



Remember also that for the Kaluza-Klein case nothing depends on the "fifth" dimension, 
and as such F^^j = 0. We get the usual Kaluza-Klein expressions, 

R^l ^ A^R^^ + ^9"e2*F^« + e2^V"F^„ = R^^, 

Rf^u ^ Rf^u + 2A(^i?^)5 - A^A^R55 - ^e^^F^^F^^ - V^O^^ - d^c^d^cj), 
R^R- ^e^^F^'^Fap - 2V"da(l) - 29„09°</). 
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6.1.2.2 S''' 



A more interesting case is performing the dimensional reduction on the 5" sphere, n = 
D — d>2. For such an isometry, the KiUing vectors ^a, a = 1, . . . , (n + l)n/2 satisfy 



(6.1.14) 



where €ab'^ are the structure constants of SO{n + 1). Because the fibre has the minimal 
dimension necessary to accommodate n(n + l)/2 independent Killing vector fields, we may 
assume without loss of generality that the Killing vector fields have components exclusively 
along the fibre: = Ca^- Furthermore, we may normalise the Killing vectors so that they 
only depend on the coordinates of the fibre, i.e. 9^^* = 0. 



Equation (6.1.3) gives the following conditions 



^ia9llU =0. 



(6.1.15) 
(6.1.16) 
(6.1.17) 



These expressions can be interpreted either as Lie derivatives of rank-2 tensors defined on 
the D-dimensional spacetime, or as Lie derivatives of a rank-2 tensor, a vector and a scalar, 
which are defined on S"". 



Together with (6.1.14), conditions (6.1.15)-(6.1.17) have the following implications: 



= fi^nM 



(6.1.18) 



because, from (6.1.15), gjj admits the maximal number of Killing vector fields and thus 



must be the metric on a maximally symmetric space at each x^. Due to (6.1.14) this 



space must be the 5" sphere, /i^" denotes the metric on an S"' with unit radius; 



9iJ.i' — 9nv{,x^) 



(6.1.19) 



because the Killing vector fields act transitively on the fibre and therefore the base 
space metric must be independent of the fibre coordinates; 



-B! = . 



because equation (6.1.16) is equivalent to 

Ka,i3^] = , 



(6.1.20) 



(6.1.21) 



and there exist no non-trivial vector fields on 5" for n > 2 that commute with all 
Killing vector fields on the sphere. 
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We write the metric on the sphere as 



(6.1.22) 



with = (j){x'^). Our D-dimensional metric has a block diagonal form. Making g^j = e^'^'h 



1] 



and = in the expressions (6.1.9)-(6.1.12) we get 



Rpiv = Rfiu - nVi^d^cp - nd^(t)dy(j), 

R = R + R- 2nV^d^(l) - n(n + \)d^ (\)d ^,(\> , 



(6.1.23) 



where R-n and R are the Ricci tensor and Ricci scalar for the metric (6.1.22). They evaluate 



to 



Rri = (n — l)/i^f, R = n{n — l)e 



-20 



(6.1.24) 



For Z)-dimensional vacuum spacetimes Rab = = = R-fy Using also (6.1.24) on (6.1.23) 
we get two coupled equations, 



e^^ {nd''(l)da(l) + V^^do^ct)) = n - 1 
These equations can also be obtained from the following action 



(6.1.25) 



1 



167rG^ 



jd / 

d x\/—ge 



R + n{n- l)e-2'^ + n(n - l)d^(t)d^'(t) 



(6.1.26) 



Performing the substitution e^*^ = will be useful for the upcoming numerical implementa- 
tion. We get 



np 



1 X-'d'^XdaX + V^daX 



2{n-l), 



R, 



(6.1.27) 



For completeness, we write in appendix |6.B the equations of motion obtained when we write 



the action (6.1.26) in the Einstein frame. 



6.2 Dimensional reduction on a (D — 4)-sphere and 3+1 split 

In the previous section we were considering a dimensional reduction under the full isometry 
group of the higher-dimensional spacetime. In the case of head-on black hole collisions, this 
would produce a reduction down to 3 spacetime dimensions. In practice, we are actually 
interested in performing a 4 -|- (L) — 4) split of the D dimensional spacetime. This may be 
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done as follows. The metric on a unit ^ may always be written in terms of the line 
element on a unit S^~^, denoted by as follows, 



hf-" 'dx'dx'^ = d9^ + sin^ 9dnD-4 ■ 



(6.2.1) 



where is a polar- like coordinate, 6 G [0, vr]. Now we introduce four dimensional coordinates, 
^ti = (^^p.^ 0-^^ /i = 0, 1, 2, 3, and define a four dimensional metric 



g^.dx^'dx'' = gppdx'^dx^ + /{x^dO^ 
as well as a new conformal factor 

X{xf') = sin^ Ogee ■ 



(6.2.2) 



(6.2.3) 



As we have seen in the previous sections, the most general Z?-dimensional metric compatible 
with SO{D — 2) isometry is, for D > 5 



ds^ = ^^^dx'^dx'' + Xixf')dnD-4: ■ 



(6.2.4) 



The geometry (6.2.4) has a manifest SO{D — 3) symmetry. We will now perform a dimen- 
sional reduction on a (D — 4)-sphere, which yields, from the D-dimensional vacuum Einstein 
equations, a set of 3+1 dimensional Einstein equations coupled to quasi-matter. In cases with 
larger symmetry (if S0{D—2) is the full isometry group, for example), the quasi-matter terms 
do not contain independent degrees of freedom and could in principle be fully determined by 
the 3 + 1 dimensional geometry. For such cases we could perform the dimensional reduction 
on a (L> — 3)-sphere instead (which has the full isometry group SO{D — 2)), which would 
yield a 2 + 1 dimensional system. The former method allows, however, the use of existing 
numerical codes, with small changes, which justifies our choice. 

The SO{D — 3) isometry group allows the study of a large class of black hole collisions with 
impact parameter and with spin: the collisions in which the two black holes always move 
on the same 2-plane and the only non trivial components of the spin 2-form are on that 
same 2-plane — see figure 6.1 With our framework we are able, therefore, to describe not 



only head-on collisions of spinless black holes but also a class of collisions for spinning black 



holes with impact parameter. As follows from the discussion of (6.1.20), the ansatz (6.2.4) 
describes general spacetimes with SO{D — 3) isometry in D > 6. We remark that the models 
with D > Q are actually the most interesting for phenomenological studies of large extra 
dimensions models (see for instance 



Taking (6.2.4) as an ansatz, we see from (6.1.26) that the D-dimensional Einstein-Hilbert 
action takes the form (for reasons related with the numerical implementation, we now use 
the variable A instead of the previously used (j)) 



1 



167rG4 



d^x^ 



g-4 

-5A 2 



R + {D-A){D-^)[x-^ + ^X-'^d^Xd^'>}j , (6.2.5) 
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Head-on 
Spinless 



SO{D - 2) 
isometry 



Non head-on 
Spinning 



SO{D - 3) 
isometry 



-^1 ^2 



^D-3 



Figure 6.1: Z?-dimensional representation, using coordinates {t,x^,x'^ 



), of two 

types of black hole collisions: (left panel) head-on for spinless black holes, for which 
the isometry group is SO{D — 2); (right panel) non head-on, with motion on a single 
2-plane, for black holes spinning in that same plane only, for which the isometry group 
is SO{D — 3). The figures make manifest the isometry group in both cases. 



where the Z)-dimensional Newton's constant G/j is related to the four dimensional one G4 
by the area of the unit D — A dimensional sphere: G4 = Gd/A^° ExpUcitly, the D- 
dimensional Einstein's equations in vacuum yield the following system of four dimensional 
equations coupled to a scalar field: 



R 



D-A 
2A 



V^5^A = 2{D - 5) - ^^d^Xdi^X . 



(6.2.6) 
(6.2.7) 



In these equations, all operators are covariant with respect to the four dimensional metric 



g^^. These could also be obtained from equations (6.1.27) withp = 1. The energy momentum 
tensor is 



D 



IQirX 



1 D — 5 

^i^dyX - —d^XdyX - {D - 5)g^u + g^ydaXd°'X 



(6.2.^ 



With this four dimensional perspective, the usual 3-1-1 split of spacetime can be performed, 



as outlined in section 2.5 As explained therein, the projection operator 7^1, and the normal 
to the three dimensional hypersurface E, {n^n^ = —1), are introduced 



(6.2.9) 



as well as the lapse a and shift j3^, 

dt = an + f3, (6.2.10) 
where t is the time coordinate. The four dimensional metric is then written in the form 



ds^ = g^udx^'dx" = -a^dt^ + 7ij(dx^ + /3Mi)(dx^' + /3Mt) , j = 1, 2, 3 . (6.2.11) 
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As usual, we introduce the extrinsic curvature Kij = —^Cnjij, which gives the evolution 

(6.2.12) 



equation for the 3-metric (2.6.2). Defining the variable 



we further get an evolution equation for A 

{dt-Cp)X = -2aKx. (6.2.13) 

Using the relation 

Da,D^X = -K^prfd^X + raYp^udf^X , (6.2.14) 
where Da denotes now the covariant derivative with respect to the 3-metric 7^j^ on S, and 



equation (|6.2.7|) we can get an evolution equation for Kx. The contraction of equation (6.2.14) 

□A = Y^D.djX - 2KKx - n^rfV^d^X . (6.2.15) 



with (7"^, yields 



Noting that 



and 



we obtain 



n^'Vau" = -D^a. 
a 



a 



Noticing also that D^adyX = ^''diadjX, we write 



□A = -f'WidjX - 2KKx + 2CnKx + -f^diadjX . 

a 

Moreover, from equation 

Df,X = 7%9,A = d^X - 2n^Kx , 

we get 

daXd'^X = -f'^diXdjX - AKI , 
so that the evolution equation for Kx is 



(6.2.16) 
(6.2.17) 
(6.2.18) 

(6.2.19) 
(6.2.20) 
(6.2.21) 



1 



{dt - Cp) Kx = -:^f^diXdja + {D - 5) + KKx + ^^Kf - ^^f^diXdjX 



2a 
1 



X 



4A 



D'^dkX. 



(6.2.22) 



Equations (6.2.13) and (6.2.22) are the evolution equations for the quasi-matter degrees of 
freedom. 
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6.2.1 BSSN formulation 



For numerical implementation, we write the evolution equations in the BSSN formulation, 
as introduced in section |4j The evolution equations are those of (4.1.14) with source terms 
determined by ( 4.1.16[ ) where the energy momentum tensor is given by equation (6.2.8). A 
straightforward computation shows that 



i7r{E + S) 
D-A 



L»-4 



D-6 



X-'^xf^diXdjX - X-^KKx - (D - 5)X-^K 



X 1 



^x^'^DidjX + {diXdjx + djXdix - I'^'iiAXda 



(6.2.23a) 



-XX "^diXdjX 



X-^KxA^ 



ll^A~'x'/^''D, [x-'/'d,x) + ^^.,X-\^'^%Xd,X, 

(6.2.23b) 



2X-^f^djKx - X-'^Kxf^djX - f^^^^AuiX-^djX - X^KX-^djX 



3.2.23c) 



where Di is the covariant derivative with respect to 7ij. 
Finally, the evolution equations for A and Kx are 

{dt -Cp)X = -2aKx, 



{dt-C^)Kx = a\{D-5) + 



6-D 



[X'^f^diXdjX - AX-^K 



+ KKx 



(6.2.24a) 
(6.2.24b) 



As stated before, in the case of head-on collisions of spinless black holes the full symmetry 



of the D-dimensional system we want to consider makes equations (6.2.24) redundant, by 



virtue of (6.2.3). This allows to determine the quasi-matter degree of freedom in terms of 



the three dimensional spatial geometry, at each time slice. The extra symmetry manifests 
itself in the fact that 7jj possesses, at all times, (at least) one Killing vector field. If one 
chooses coordinates adapted to this Killing vector field, d/dO^ the metric can then be written 



in the form (6.2.2), and then the quasi-matter degree of freedom can be determined from the 



spatial geometry by (6.2.3). In the numerical implementation, one can either determine, at 



each time-step, the scalar field through (6.2.3), or impose (6.2.3) only in the initial data, and 



then evolve the scalar field using equation (6.2.24). We have implemented the latter method 



6.3 Higher-dimensional initial data 

Having written our system of evolution equations, we now need to construct relevant initial 
data. Building on the results outlined in chapter [3] (based on |1271 1128j ) , we now present a 
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generalisation of the spectral solver in [139] that generates initial data for black hole binaries 
in D > 5 dimensions with non- vanishing initial boost [177j . 

In this section we need to make some changes to our notation. Early lower case Latin indices 
a, b, c, . . . will here extend from 1 to I? — 1, late lower case Latin indices i, j, k, . . . run from 
1 to 3 and early upper case Latin indices A, B, C, . . . from 4 to D — 1. 



6.3.1 Coordinate transformation 

We start by recalling, from chapter [3] that, for a system of boosted black holes, we can 



solve the momentum constraint equation (2.6.3d) analytically. It remains then to solve 



equation (3.2.35), which we re- write here: 



The numerical solution of this equation will be our task in this section. 

First, it is convenient to transform to a coordinate system adapted to the generalised axial 
symmetry SO{D — 2) in D = 5 dimensions and SO{D — 3) in D > 6 dimensions as discussed 
in section |6.2[ For this purpose we consider the (flat) conformal spatial metric in cylindrical 
coordinates 

^ah<^x''dx^ = dz'^ + dp^ + {dip"^ + sin^ipdinD~i) , (6.3.2) 

where is the metric on the [D — 4)-sphere. Observe that 99 is a polar rather than an 

azimuthal coordinate, i.e. G [0,7r]. Next, we introduce "incomplete" Cartesian coordinates 
as 

x = pcos(/3, y = psinif ^ (6.3.3) 

where —00 < x < +00 and < y < +00. The D dimensional initial data for the spatial 
metric is then 

7afcdx"dx'' = [dx^ + dy2 + dz^ + y'^dVto-i] ■ (6.3.4) 

We can transform the D — \ dimensional Cartesian coordinates X"" = (x^, . . . ,x^~^) to the 
coordinate system = (x, y, ^1, ^2, • • • , io-A) with hyperspherical coordinates ^1, . . . , ^d-a 

by 

X^ = X 

x^ = y cos^i 
X = z 



= y sin^i cos^2 {D>Q) 

y sin ^1 sin ^2 cos ^3 (£> > 7) . (6.3.5) 



x^ 



x^ ^ = y sin ^1 • • • sin cos iD>7) 
x^~^ = y sin ^1 • • • sin cos {D > 6) 

x^~^ = y sin ^1 • • • sin {D > 5) 
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Without loss of generality, we can always choose coordinates such that the black holes are 
initially located on the z axis at zi and Z2 and have momenta of equal magnitude in opposite 



directions P^-^^ = ~P(2)- Inserting the momenta into equation (3.2.23) then provides the 
conformal traceless extrinsic cuvature and the differential equation (6.3.1) which is solved 
numerically for u. 

The class of symmetries covered by the formalism developed in this chapter includes head-on 
and grazing collisions of non-spinning black holes with initial position and momenta 

x^,) = (0,0,zi,0,...,0), x^2) = (0,0,Z2,0,...,0) 

= (P^ 0, P^ 0, . . . , 0) = -Pf2) . (6.3.6) 

Note that a non-zero P^ is not compatible with the assumed symmetries. On the other hand, 
the X-axis can always be oriented such that the collision takes place in the xz plane. Our 
formalism therefore covers general grazing collisions of non-spinning black hole binaries in D 
dimensions. 



6.3.2 Four dimensional initial data for a general D head-on collision 

We will now discuss in detail the case of black holes with momenta in the z direction, that 



is, the case given by setting P^ = in equation (6.3.6). The linear momenta are thus given 

by 

= (0, 0, P^ 0, . . . , 0) = -P^^y (6.3.7) 



The rescaled trace-free part of the extrinsic curvature for such a configuration is 

Aab 



^ab + ^ab ' 



(6.3.8) 



where A^^^ and A^^^ are given by equation (3.2.23) with (6.3.6) and (6.3.7). Using equa- 



tion (|6.3.5|) we can write this in the coordinate system y"" adapted to the spacetime symmetry: 

(6.3.9) 



47r(L> - 1)P^ / 







{D - 2)Ad-2{x^ + y^ + {z- zi)2)^ ' 


\ 


"■AB 



with 

(1) _ 



a, 



and 



-^ab 



-(D~4:)x'^+y^+(z-zxf](z-zx) 
{D-3)xy{z-zi) 
a;[z2+y2 + (D-2)(2-2i)2] 



{D-3)xy{z-zi) a;[x2+y2+(D_2)(z-2i)2] \ 

-[x^-(D-4)y^+(z-zi)^]{z~zi) y[x^+y^+(D-2)(z-zi)^] 

S/[x2+y2 + (D-2)(2-2i)2] [xHy^ + iD-2){z-zi)^]iz~-zi) J 

(6.3.10) 



^AB = -y^i^ - ^i) [^^ + + zif] hAB 



(6.3.11) 



where h^s is the metric on the {D — 4)-sphere. The expression for is analogous, but 



with Z2 in place of z\ and —P^ in place of P^ in equation (6.3.9). 



We now need to re-express these quantities in terms of our 3 + 1 quantities, (7ij, Kjj, A, K\)^ 
as introduced in the previous section. These are the variables evolved in time and therefore 
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the variables we ultimately wish to construct from the initial data calculation. For their 
extraction we first note that jij, Kij and Kx are related to the {D — l)-dimensional metric 
7aft and extrinsic curvature Kab by 

(6.3.12) 

, Kab = -^KxhAB , 
D -AK\ 

K = K + — ^. (6.3.13) 





= lij 


liA 


= , 






KiA 


= , 



Using these relations and equation (6.2.4) we can express all "3+1" variables in terms of 
those describing the initial data 



K,, = ^-2(ilj^ + igO , Kx = 2i;-V(.P^ + P-), (6.3.14) 
{D-4)Kx 



where 



2A 

4tt{D - l)P%z - zi] 



(D - 2)Ad-2{x^ + + {z - 
p- ^ ATr{D-l)P'{z-Z2) 

' (D- 2)AD-2ix^ + 2/2 + (z - Z2)2) V 

The conformal factor is 

_Mi ^ /f2 

4 [x2 +y2 + {z- Zi)2](^-3)/2 ^ + y2 + _ ^2)2] 

and u is the solution of the equation 



0-1 ' 

I 2 



(6.3.15) 



~ ^ . [.2 , „,2 , ^. ^2^{^-3)/2 ^ . r^2 , „2 , _ .„^2l{^-3)/2 + ^ ' (6.3.16) 



D - 3 \ 3-D - - 3g-5 



where 

A'^'Aab = {A\f + ijf )(i*^'(i) + i*^'(')) + (Z) - 4)(P+ + p-f . (6.3.18) 



Our numerical construction of the function u will be based on the spectral solver developed 
in [139j . This solver employs coordinates specifically adapted to the asymptotic behaviour of 
u at spatial infinity. In order to investigate this behaviour, we next consider a single black 
hole with non-zero linear momentum. 



6.3.3 Single puncture with linear momentum 



For a single puncture with momentum P^ located at the origin z = 0, equation (3.2.23) 
implies 



A'^'^A, 



ab 



16-K^{D - 1) 



{D - 2YAl_^r^^D-2) 



pz 



2 + D{D-3) 



(6.3.19) 
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so that equation (6.3.17) takes the form 
87r2(D-l)2p-3) . 



An + 



1 + 



L>(D - 3) /z\2 



3-D-5 

D-3 = . 



(6.3.20) 



It turns out to be convenient for solving this differential equation to introduce a hyperspherical 
coordinate system on the D — 1 dimensional spatial slices, such that the flat conformal metric 
is 

ds^ = 7a6dx"dx^ = dr^ + [M^ + sin^"^ {dcp"^ + sin^ ipdQD-i)] , 



with cos'd = -. We further introduce the radial coordinate 

r 



X 



-1 



(6.3.21) 



which reduces to the coordinate A of equation (31) in [139j for the case of D = 4 spacetime 



dimensions. Expressed in the new coordinate system, equation (6.3.20) becomes 



dxx + -Tpdx + 



1 



(D- 3)2X2(1 -X)2 



d^^ + (D — 3) cot -dd^ 



+ . o ^ (i9w + {D - A) cot ipd^) }u 
sm v 

P.V ..-o^ . „._3o^ / D(D-3) 

D-'i 1 J ^ ^ 



-a I — ) X {l + uX) 



+ ^cos2i?) , (6.3.22) 



with 



a = 



1287r2(D - 1) 



{D-?,){D-2YAl 



D-2 



For L) = 4 we recover equation (40) of |139| . In order to study the behaviour of the solution 
at spatial infinity, we now perform a Taylor expansion m v = — , 



(6.3.23) 



Odd powers of v have to vanish in order to satisfy equation (6.3.22). We have the following 
equation for ui 



dxx + -T?dx + 



X^ (D - 3)2^2(1 -X) 



2 [dm + {D -3) cot -dd^] I ui 



-aX-^ (l + cos2 7?) . (6.3.24) 



In order to solve equation (6.3.24), we make the ansatz 



ui = f{X)+g{X)QD{cos^) , 



(6.3.25) 
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where QoCcost?) = — 1) cos^ ■!?- 
f{X) and g{X) take the form 



1. By solving equation (6.3.24), we find that the functions 



fix) 



327r2(D 



3) 



D+i 



(6.3.26) 



giX) = h 



X 



a 



1 -X, 

D{D-3f 
2(D+1)(D-1) 



2 D-1 

D-s fl-X\D-a 
+ k2 



1 

D-1 



1 - X)d 

^X^{1- X)'b^2Fi' 2 



-2-^1 



D-1 D-1 ■ o D-2 . y 

D-3' D-3' ^D-3' 



2D 



D-3' D-3' '^D-3' 



(6.3.27) 



where 2Fi{a,b; c; X) is the hypergeometric function and ki^2 are constants to be fixed by 
imposing that g{X = 1) = and g{X = 0) is smooth. Requiring analyticity at X = and 
using the property F{a, b, c, 0) = 1, we immediately find k2 = 0. 

We are now interested in the large X — )• 1 limit. Therefore, we use the z ^ 1 — z 
transformation law for the hypergeometric functions |186| . 



F{a-c+l,b-c+l,2-c, z) = 

r(2-c)r(a + 6-c) 



{1-zY 



+ 



r(a-c + l)r(6-c + l) 
r{2 -c)T{c- a -b) 



F{l-a,l-b,c-a-b+l,l-z) 



r(l-a)r(l-6) 
Requiring a regular solution we find that ki has to satisfy 



F{a-c+l,b-c+l,-c+a+b+l,l-z) . (6.3.28) 



k. 



64tt^D{D - 3) 



2 r 



2(D-2) y 
D-3 ) 



{D-2nD + l)Al_, rf3D 



D-3 



(6.3.29) 



Let us write these functions explicitly for D = 4, 5, 7 (for D = 6 the hypergeometric function 
does not simplify): 



D = 4 : 



f{X) = -{l-X') , 



(6.3.30) 



10X3 



84(1 - X) log(l -X)+ 84X - 42X2 _ -^^^3 _ _ _ 

(6.3.31) 



These are equations (42-44) in |139| . with appropriate redefinitions. 



CHAPTER 6. HIGHER-DIMENSIONAL NUMERICAL RELATIVITY 76 



D 



D 



fix) 

9iX) 



9iX) 



128 d - X^l 



80(1 -X 
8l7r2X2 



l2 r 



4 log(l - X) + 4X + + 



(6.3.32) 
(6.3.33) 



257r4 



28 



1257r4^(l - X)X^ 

+ SttX^ + 6 (5 - lOX + 4X2) ^^^g-j^ ^ 



(6.3.34) 

30V(1 -X)X + 40y^{l-X)X3 - 16\/(1 -X)X7 

(6.3.35) 



Analysing these expressions, we can anticipate the convergence properties of the numerical 
solutions obtained in terms of pseudo-spectral methods. For instance, analyticity of / and g 
suggests exponential convergence. As will become clear in the next section, we are interested 
in the convergence properties in a coordinate A behaving as A ~ ^ ~ ^'^^ large r. We thus 
introduce a coordinate A that satisfies 



X = (1 + {A-^ - 1)^-^) 



(6.3.36) 



In terms of the A coordinate, we find that the functions / are analytical. For the function g 
in the vicinity of j4 = 1, the leading terms behave as follows: 



D = 5 



D = 6 



D = 7 



9{A) 



80 
81^ 



(l-A)^[81og(l-^) + 7] , 



19683 ,^ ,,5 
6272^(1-^) 



84 
257r3 



(1 - Af 



(6.3.37) 



(6.3.38) 



(6.3.39) 



From the behaviour of the functions / and g and equation (6.3.25) we conclude that the 
first term in the expansion (6.3.23) has a leading-order behaviour ui ~ l/r^~^ 



as r 



00. 



Iteratively solving equation (6.3.22) for higher powers of v is complicated by the presence of 



the source terms on the right hand side, but under simplifying assumptions indicates that 
higher-order terms Uj > 2 acquire additional factors of 1/r and therefore the leading-order 
fall off behaviour is given correctly by that of ui . This result is confirmed by our numerical 
investigation using finite boost parameters as we shall discuss in the next section. 
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With regard to the analyticity of the solutions and the resulting expectations for the con- 
vergence properties of a spectral algorithm, we summarise the results of our analytical study 
of a single puncture as follows. In D = 6, 7, the leading terms are analytic functions in the 
vicinity of A = 1. Actually, for D = 7, g{A) is analytic in the vicinity of any point. Therefore, 
we expect exponential convergence of the pseudo-spectral code. For D = 5, one observes the 
presence of a logarithmic term. This type of term is known to arise in Z) = 4, when punctures 
have non- vanishing momenta \187\ I188| and in that case their presence makes the convergence 
algebraic in the single puncture case. In the next section we shall investigate the impact of 
the logarithmic terms on the convergence properties of our spectral solver. 



6.3.4 Two punctures with linear momentum 



6.3.4.1 Code changes 



We first explicitly list the modifications applied to the spectral solver of reference |139j and 
demonstrate how these modifications enable us to generate initial data for boosted black 
hole binaries with convergence properties and levels of constraint violation similar to the 
-D = 4 case. For this purpose we start by recalling that the spectral solver of [139| employs 
coordinates 

^G[0,1], Be [-1,1], 0G[O,27r], (6.3.40) 
which are defined by equation (62) of |139j . 

, 2A 1-B^ 

42 + 1 2B 



where b is half of the coordinate distance between the punctures. In particular, the coordinate 
A satisfies 

r ^ 00 A^ I . (6.3.42) 
The first modification consist in adapting the source term and Laplace operator according 



to (6.3.17) 



Next, we note that the type of high-energy collisions which form the main motivation for 
this work often start from relatively large initial separations of the holes, l^i — 2;2| ^ fs- In 
order to obtain high-precision solutions for such binary configurations, we found it crucial to 
introduce a coordinate A' defined as 

, sinhk(A' + l)/2l 

A = ^ \ ' , (6.3.43) 

smh K 

where k is an adjustable free parameter. Note that for k = we obtain A = ^{A' + 1) For 
K > 0, however, the new coordinate A' provides the spectral method with enhanced resolution 
near ^ ~ 0. 
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Table 6.1: ADM mass obtained with equation (6.3.47) in units of the "bare" Schwarzschild radius 



„D-3 



„D-3 



„D-3 



+ The variation of the ADM mass with resolution is of the order of 

10^ for all D and n > 100 grid points indicating that the accuracy in the ADM mass 
is limited by round-off errors. 



D 



b/rs 



P/^ 



,D-3 



,D-3 



MADM/rg' 

1.78 
2.27 
2.96 
3.81 



4 
5 
6 
7 



30.185 
30.185 
30.185 
30.185 



0.8 
0.8 
0.8 
0.8 



3.555 
1.931 
1.415 
1.236 



A further modification is related to the asymptotic fall off of the function u as obtained in 
the previous section, 



1 



u 



rD-3 



(6.3.44) 



To naturally accommodate this behaviour with the spectral coordinates used in the code, we 
have changed the variable U of equation (5) in [139] to 



u = {A'- i f-'^U . 
Note that this U variable is the variable that the code actually solves for. 



(6.3.45) 



Finally, we adjust the calculation of the ADM mass from the numerical solution. For this 
purpose, we note that, asymptotically 



Ip = 1+ n Q + n Q + U ~ 1 + 



4r^-3 4rf^-3 



(6.3.46) 



with a = ro~^ = is^^adm ^ = rj.^^ ^j^jyj ^^^^ -g ^^len obtained from 

^ -^global AD-2(iJ—^) ^ •'(±) 



r , = rc, ^ + r-o ^ + 4 lim 



'-'global "Ji 



(+) 



-2b- 



'^f^y'\{A'=i), 



(6.3.47) 



where we have used equation (62) of |139]. and equation (6.3.43) and (6.3.45). We show in 
table [nH] the values obtained for the ADM mass of some cases we considered. 



6.3.4.2 Results 

We now study the numerical results as obtained for D = 4,5,6,7 with these adaptations of 
the spectral solver of |139j . Throughout the remainder of this section we will graphically 
present results in units of the "bare" Schwarzschild radius defined as rc.~ =?'c~ -\-ra~ ■ 

We first address the convergence properties of the numerical algorithm by evaluating the 
quantity 

Sn,miu) = max|l - Un/Um\ , (6.3.48) 
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where the maximum is obtained along the colhsion axis, i.e. z-axis in our case. Here, the 
index m refers to a reference solution obtained using a large number m of grid points 
while n denotes test solutions using a coarser resolution, n < ni. The result obtained 
for black hole binaries with initial separation b/rs = 30.185 and boost P^/rg~^ = 0.8 in 

= 4, 5, 6 and 7 dimensions is displayed in figure 6.2 
achieving a given target accuracy 6n,r, 



We note from this figure, that 



requires a larger number of points n as D increases. 
We emphasise in this context, however, that this increase in computational cost in higher 
dimensions is unlikely to significantly affect the total computational cost of the simulations 
which typically are dominated by the time evolution rather than the initial data calculation. 
Most importantly, we observe exponential convergence up to a level of 6n,m (u) « IQ-^ for 
all values of the spacetime dimensionality D. Below that level, the two leftmost curves in 
figure 6^ corresponding to D = A and D = 5, respectively, show that the rate of convergence 
decreases indicating that the logarithmic terms become significant and reduce the convergence 
to algebraic level similar to the observation in figure 4 of reference |139j . For D = 6, the 
convergence remains exponential, in agreement with the absence of logarithmic terms in 
the analysis of section [6. 3. 3[ Irrespective of a change to algebraic convergence, however, our 
algorithm is capable of reducing the quantity 5m,n {u) for all values of -D to a level comparable 
to the case D = 4 and, thus, producing initial data of similar quality as in 3+1 dimensions, 
provided we use a sufficiently high resolution n. 
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Figure 6.2: Convergence plot for the b/rs = 30.185, P/r§-^ = ±0.80 cases. 



For illustration, we plot in figure 6.3 the function u obtained for the case of b/rs = 30.185, 
P^/r^~^ = 0.8. The behaviour is qualitatively similar for all values of D, but the figure 
demonstrates the faster fall off for larger D as predicted by (6.3.44). For this plot we have 
used = 300, tlb = 300 and = 4 grid points. The inset in the figure shows the function 
u in the immediate vicinity of the puncture. While the profile develops multiple extrema for 
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D > 4:, the profile remains smooth for all values of D. 



10' 
10" 
10"' 
10"- 
10" 
10"' 
10"' 
10" 
10"' 
10' 

10 

10 



P/rr 



= 0.8, b/rs =30.185, n=300 




10 




Figure 6.3: u function ioi D — A, . . . ,7 plotted along the 2:-axis, in units of rs- We used ua = ns 
n — 300, rij, — 4. We also show a zoom around the puncture. 



Finally, "we sho"w in figure 6.4, the Hamiltonian constraint corresponding to the solutions 
presented in figure |6.3| as measured by a fourth-order finite differencing scheme of the Lean 



evolution code (section 4.3 and reference |149j ). We emphasise that the violation of equa- 
tion (6.3.1 ) inside the spectral initial data solver is < 10~^^ by construction. The independent 



evaluation of the constraint violation in the evolution code serves two purposes. First, it 



checks that the differential equation (6.3.1 ) solved by the spectral method corresponds to the 



Hamiltonian constraint formulated in ADM variables; an error in coding up the differential 



equation (6.3.1) could still result in a solution for u of the spectral solver, but would manifest 



itself in significantly larger violations in figure 6.4 Second, it demonstrates that the remaining 
numerical error is dominated by the time evolution instead of the initial solver. Note in 
this context that the relatively large violations of order unity near the puncture location in 



figure 6.4 are an artifact of the fourth-order discretisation in the diagnostics of the evolution 
code and are typical for evolutions of the moving-puncture type; see e.g. the right panel in 
figure 8 in Bro'wn et al. [189]. 

The solid (blue) curve obtained for the "standard" = 4 case serves as reference. For all 
values of D the constraint violations are maximal at the puncture location zi/rs ~ 15 and 
rapidly decrease away from the puncture. As expected from the higher fall off rate of the 
grid functions for larger D, the constraints also drop faster for higher dimensionality of the 
spacetime. 
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= 0.8, b/rg =30.185, n=300 




Figure 6.4: Violation of the Hamiltonian constraint along the z-axis, evaluated with a fourth order 
finite difference scheme. The growth of the constraint violation near the puncture is an 
artifact of finite-differencing across the puncture; see text for details. 

6.4 Numerical evolutions 



Having established that our initial data code is working, we will now show some numerical 



results, obtained by adapting the Lean code introduced in section 4.3 In this section we 
will begin by briefly commenting on numerical issues generated by the quasi-matter terms 
arising from the dimensional reduction. We then present some code tests and results. 



From the initial data construction of section [6^ we see that the quasi-matter field A has a 
fall off as y — )• 0, that is, on the xz plane (cf. (6.3.14)). From (6.2.23), we see that this 
leads to divisions by zero on the right-hand side of the BSSN evolution equations; thus, we 
need to isolate such irregular terms and re-write the equations in terms of variables which 
are explicitly regular at y = 0. In this spirit, we introduce the following evolution variable 

X 



and corresponding auxiliary variable 
1 

2ay^ 



K, 



C 



1 

'2a 



(6.4.1) 



(6.4.2) 



The full quasi-matter terms and evolution equations in terms of these regular variables can 
be found in Appendices A and B of |175j . Here we note only that, with the above definition, 
we have the relation 

1 



y 



X 3 X 



(6.4.3) 
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Finally, for long term evolutions, we employed the following gauge conditions, which are 



generalisations of conditions (4.2.3), (4.2.6) 



dt - P'dk) a = -2a{r,KK + r,K,K^) , 



(6.4.4) 
(6.4.5) 



Note the extra term involving in the slicing condition compared with standard moving 
puncture gauge in 3 + 1 dimensions and the additional freedom we have introduced in the 
form of the parameters tjk and r]K^ ■ 



6.4.1 Code tests 
6.4.1.1 Geodesic slicing 

As a first test of our numerical implementation, we have numerically evolved a single D = 5 
Tangherlini black hole in the so-called geodesic slicing, which corresponds to fixing the gauge 



parameters to (4.2.1) throughout the evolution. Such a gauge choice is not adequate to 
perform long term numerical evolutions. The advantage of this choice, though, is that one 
can easily write the Tangherlini metric element in this coordinate system, which we can 
then match against the numerically obtained solutions [971 . This coordinate system may 
be achieved by setting a congruence of in- falling radial time-like geodesies, each geodesic 
starting from rest at radial coordinate rg, with rg spanning the interval [/i, +oo[, and using 
their proper time r and tq as coordinates (instead of the standard t, r Schwarzschild-like 
coordinates). The line element becomes 



2 



where rQ{R) is given by 

ro(i?) = r{i + ^Y (6-4-7) 

Before the breaking down of the numerical evolution, we can compare our numerical results 
with the above metric element. This is shown in figure [63} where we have plotted one metric 
component 'yxx along the x axis (left) and C,/x (right), for various values of r using both the 
analytical solution and numerical data. The agreement is excellent for and good for C/x- 
The latter shows some deviations very close to the puncture, but we believe that it is not a 
problem for two reasons: (i) the agreement improves for higher resolution; (ii) the mismatch 
does not propagate outside of the horizon. 
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x//i y//i 



Figure 6.5: Numerical values versus analytical plot oi '^xx along the a;-axis (left panel) and of C/x = 
X/y^ along the y axis (right panel), for various values of r, for the single Tangherlini 
black hole in five dimensions. The horizontal axis are in units of /i. 



6.4.1.2 Single black hole evolution 



To further test our numerical framework, we have performed long term simulations of a single 
black hole m D = 5 using the gauge conditions in (6.4.4) and (6.4.5), the initial data from 



equations (6.3.14) (with P"*" = = P ) and grid setup (cf. section 4.3) 



{(512, 256, 128, 64, 32, 16, 8, 4, 2) x (), h} , 



in units of fx with resolutions he = 1/32 and /if = 1/48. In figure 6.6 we show the Hamiltonian 
constraint and the y component of the momentum constraint at evolution time t = 28/U. 
For the Hamiltonian constraint the convergence is essentially 4th order; for the momentum 
constraint it decreases towards 2nd or 3rd order in patches. 




y/fi 



y/li 



Figure 6.6: Hamiltonian constraint (left panel) and y-component of the momentum constraint (right 
panel) at time t — 28/i, for the evolution of a single Tangherlini black hole in five 
dimensions. 
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6.4.1.3 Head-on collision 



Finally, we tested the code capability to evolve a head-on collision from rest. Using again the 

= , we let two black holes with parameters 

(6.4.8) 



initial conditions (6.3.14) with P 



/f2 

2 ' 

zb = 3.185 /i . 



).4.c 



collide from rest, using the grid setup 

{(512, 256, 128, 64, 32, 16, 8) x (2, 1), h = 1/32} , 



in units of fi. The gauge variables a and /?* were evolved according to equations (6.4.4) and 
(6.4.5) with parameters rjK = VKi; = 1-5 and r] = 0.75. 
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Figure 6.7: The BSSN variable x (l^ft panel) and the quasi-matter momentum K^^ (right panel) 
shown along the axis of collision for a head-on collision at times i = 0, 5, 20, 40 and 
256 ^. Note that iiT^; = at i = 0. 



In figure 6.7 we show the conformal factor x and the momentum along the axis of collision 
at various times for such an evolution. At early times, the evolution is dominated by the 
adjustment of the gauge (cf. the solid and short-dashed curves). The two holes next start 
approaching each other (long-dashed and dotted curves) and eventually merge and settle down 
into a single stationary hole (dash-dotted curves). No signs of instabilities were observed. 



6.4.2 Head-on collisions 

Having established, in the previous section, that our numerical implementation does work 
(for the five dimensional case) we will now present some results. 

Since the final result of a head-on collision of two D dimensional, non-spinning black holes 
approaches, at late times, a D dimensional Schwarzschild (i.e. Tangherlini) black hole, we 
can make use of the Kodama-Ishibashi formalism, presented in section |5.1.2| to extract the 
gravitational wave information. Our remaining task is then to obtain the relevant gauge- 
invariant quantities from our numerical data. We will give here the main steps for this 
procedure; the full details can be found in reference |176j . 
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6.4.2.1 Coordinate frames 



In the approach developed in section 6.2 we perform a dimensional reduction by isometry on 
the {D — 4)-sphere in such a way that the D dimensional vacuum Einstein equations 

are rewritten as an effective 3 + 1 dimensional time evolution problem with source terms 
that involve a scalar field. We focus here on Z? > 5 dimensional spacetimes with SO{D — 2) 
isometry group, which allows us to model head-on collisions of non-spinning black holes; we 
dub hereafter these spacetimes as axially symmetric. Although the corresponding symmetry 
manifold is the (D — 3)-sphere S^~^, the quotient manifold in our dimensional reduction is its 
submanifold S^~^. The coordinate frame in which the numerical simulations are performed 
is 

(x^ ct>\..., c^""-^) = (t, x,y,z,^\..., c/>^-^) , (6.4.10) 

where the angles (p^ , . . . , (f)^~'^ describe the quotient manifold S^~^ and do not appear 
explicitly in the simulations. Here, z is the symmetry axis, i.e. the collision line. 



Recall that in the frame (6.4.10), the spacetime metric has the form (cf. equation ( 6.2.4[ )) 



ds^ = -a^dt^ + -fij{dx' + /3Mt)(dx^' + (3^dt) + X{x^')dnD-4 , 



(6.4.11) 



Also recall, from equations (6.2.2) and (6.2.3), that with an appropriate transformation of the 



four dimensional coordinates x^, the residual symmetry left after the dimensional reduction 
on S^~^ can be made manifest: x^^ — )■ (x'^, 9) (// = 0, 1, 2), 



5^,(x")dx'^dx'^ = 5^p(x")dx'^dx^ + g0e{x'')de^ 



(6.4.12) 



and 



A(x^) = sin2 Ogeeix^) 



(6.4.13) 



so that equation (6.4.11) takes the form ds = g^^dx^^dx'^ + gg^dQc-s- 



To extract the gravitational waves with the KI formalism, spacetime, away from the black 
holes, is required to be approximately spherically symmetric. In D dimensions this means 
symmetry with respect to rotations on S^~'^, which is manifest in the coordinate frame: 



X" 



-D-4 



(6.4.14) 



Note that x"" = t,r and that we have introduced polar-like coordinates 6,6 £ [0, vr] to "build 
up" the manifold S^^"^ in the background, together with a radial spherical coordinate r, 
which is the areal coordinate in the background. 



The coordinate frame (6.4.14) is defined in such a way that the metric can be expressed as 
a stationary background (ds*-"-*)^ (i-e., the Tangherlini metric) plus a perturbation (ds^^^)"^ 



which decays faster than ^ for large r, and the formalism from section 5.1.2 can thus 

be applied [176]. 
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6.4.2.2 Implementation of axisymmetry 

In an axially symmetric spacetime, tiie metric perturbations are symmetric with respect 
to S^^^. Therefore, the harmonics in the expansion of Hmn depend only on the angle 9 
(which does not belong to S^~^). Furthermore, since there are no off-diagonal terms in the 
metric, the only non-vanishing g^j components are g^g] the only components g^j are either 
proportional to Jij, or all vanishing but g^g. This implies that only scalar spherical harmonics 
can appear in the expansion of the metric perturbations. Indeed, if 



(V^0,...,0), V*=V*(^) 



then equation (5.1.35) gives 



Y'j = V|^^ = ^ V'^ = 
Similarly, from equation (5.1.32) we obtain T^;; = 0. 



V = 0. 



(6.4.15) 
(6.4.16) 



The scalar harmonics, solutions of equation (5.1.38) and which depend only on the coordinate 

as discussed in references |17H 1190^ 



9, are given by the Gegenbauer polynomials C^^ '^)/'^^ 
I127j : writing explicitly the index /, they take the form 

where the normalization K^^ is chosen such that 



V2cp-3)/2(cos0-), 



(6.4.17) 



(6.4.18) 



and A;2 = /(/ + £)_ 3). 

Metric perturbations, and corresponding gauge-invariant functions, can then be computed in 
terms of these functions |176j . 



6.4.2.3 Extracting gravitational waves 

In the KI framework, the emitted gravitational waves are described by the master function 

ute direc 

$i = (D-2)r(^-4)/2- 



cf. section 5.1.2 We can compute directly ^ t with |156[ 1176 

-F\ + 2rF, 



k"^ -D + 2 + 



(D-2){D-l) ' 



(6.4.19) 



where k"^ = l{l + D - 3). The energy flux can then be computed from expressions (^5.1.45^ 



(5.1.46). 



6.4.3 Head-on collision from rest in D 



Having introduced and tested our formalism and numerical code, we now present results 
obtained for head-on collisions of five-dimensional black holes. The black holes collide 
*Note that there is a factor r missing in equation (3.15) of reference |156| . 
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Table 6.2: Grid structure and initial parameters of the head-on collisions starting from rest in D = 5. 

The grid setup is given in terms of the "radii" of the individual refinement levels, in units 
of rs, as well as the resolution near the punctures h. d is the initial coordinate separation 
of the two punctures and L denotes the proper initial separation. 



Run 


Grid Setup 




d/rs 


L/rs 


HD5a 


{(256,128,64,32,16,8,4) x (0.5,0.25), h 


= rs/8A} 


1.57 


1.42 


HD5b 


{(256,128,64,32,16,8,4) x (0.5,0.25), h 


= rs/8A} 


1.99 


1.87 


HD5c 


{(256, 128, 64, 32, 16, 8, 4) x (1, 0.5), h = 


rs/M} 


2.51 


2.41 


HD5d 


{(256, 128, 64, 32, 16, 8, 4) x (1, 0.5), h = 


rs/84} 


3.17 


3.09 


HD5e 


{(256, 128, 64, 32, 16, 8) x (2, 1, 0.5), h = 


r5/84} 


6.37 


6.33 


HD5f 


{(256, 128, 64, 32, 16, 8) x (2, 1, 0.5), h = 


r5/84} 


10.37 


10.35 



from rest, with initial coordinate separation d. Note that in five spacetime dimensions the 

Schwarzschild radius is related to the ADM mass M via 

n 8M , , 

r| = — . (6.4.20) 

We therefore define the "total" Schwarzschild radius rs such that '^l = i + 2- using 
this definition, rs has physical dimension of length and provides a suitable unit for measuring 
both results and grid setup. 



As summarised in table 6.2, we consider a sequence of binaries with initial coordinate 
separation ranging from d = 3.17r5' io d = 10.37r5. The table further lists the proper 
separation L along the line of sight between the holes and the grid configurations used for 
the individual simulations. 



6.4.3.1 Newtonian collision time 



An estimate of the time at which the black holes "collide" can be obtained by considering 
a Newtonian approximation of two point particles in D = 5. The Newtonian time it takes 
for two point-masses (with Schwarzschild parameters rs^i and rs,2) to collide from rest with 
initial distance L in D dimensions is given by 



ifree-fall 2^ / L \ 2 



rs D-3 \rs 



(6.4.21) 



where rg ^ = r^^ ^ + r^ .^ ^ and 



^d. = V^£ii±4^ . (6.4.22) 

1-^ r(i + ^) ^ > 



For D = 4, one recovers the standard result tfree-faii = f y -^^/''^I'^S' > whereas for Z? = 5 we 
get 

tiree-l.ll = {L/rsfrs. (6.4.23) 
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In general relativity, black hole trajectories and merger times are intrinsically observer de- 
pendent quantities. For our comparison with Newtonian estimates we have chosen relativistic 
trajectories as viewed by observers adapted to the numerical coordinate system. While the 
lack of fundamentally gauge invariant analogues in general relativity prevents us from deriving 
rigorous conclusions, we believe such a comparison to serve the intuitive interpretation of 
results obtained within the "moving puncture" gauge. Bearing in mind these caveats, we 
plot in figure [6^ the analytical estimate of the Newtonian time of collision, together with the 
numerically computed time of formation of a common apparent horizon. Also shown in the 
figure is the time at which the separation between the individual hole's puncture trajectory 
decreases below the Schwarzschild parameter rs- The remarkable agreement provides yet 




Figure 6.8: Estimates for the time it takes for two equal-mass black holes to collide in D = 5. 

The first estimate is given by the time icAH elapsed until a single common apparent 
horizon engulfs both black holes (diamonds), the second estimate is obtained by using the 
trajectory of the black holes, i.e., the time itraj at which their separation has decreased 
below the Schwarzschild radius (circles). Finally, these numerical results are compared 



against a simple Newtonian estimate, given by equation (6.4.23) (blue solid line). 



another example of how well numerically successful gauge conditions appear to be adapted 
to the black hole kinematics. 



6.4.3.2 Waveforms 

Let us now discuss the gravitational wave signal, extracted with the KI formalism, generated 
by the head-on collision of two black holes in five dimensions. 



In figure 6.9, the / = 2 multipole of the KI function ^>^j for model HD5e obtained at different 



extraction radii is plotted. A small spurious wavepulse due to the initial data construction 
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Figure 6.9: Left panel: The / = 2 component of the KI waveform for model HD5e extracted at radii 
Rcx/rs = 20, 40 and 60 and shifted in time by Rc^/rs- Right panel: The / = 2 and I = 4 
mode of the KI function for the same simulation, extracted at Rc^/rs = 60. For clarity, 
the I — 4 component has been re-scaled by a factor of 100. 



is visible at At « 0, the so-called "junk radiation". The physical part of the waveform is 
dominated by the merger signal around At = 50rs, followed by the (exponentially damped) 
ringdown, whereas the infall of the holes before At = AOrs does not produce a significant 
amount of gravitational waves. Comparison of the waveforms extracted at different radii 
demonstrates excellent agreement, in particular for those extracted at Rex = 40rs' and GOrg. 
Extrapolation of the radiated energy to infinite extraction radius yield a relative error of 5 % 
at -Rex = 60r5, indicating that such radii are adequate for the analysis presented in this work. 

Due to symmetry, no gravitational waves are emitted in the / = 3 multipole, so that / = 4 
represents the second strongest contribution to the wave signal. As demonstrated in the 



right panel of figure |6.9[ however, its amplitude is two orders of magnitude below that of the 
quadrupole. 

In order to assess how accurately we are thus able to approximate an infall from infinity, we 
have varied the initial separation for models HD5a to HD5f as summarised in table |6.2[ As 



demonstrated in figure 6.10, for models HD5e and HD5f we can safely neglect the spurious 
radiation as well as the impact of a finite initial separation, provided we use a sufficiently large 
initial distance d > 6rs of the binary. Here, we compare the radiation emitted during the 
head-on collision of black holes starting from rest with initial separations G.STrg and lO.STrg. 
The waveforms have been shifted in time by the extraction radius i?ex = QOrs and such that 
the formation of a common apparent horizon occurs at At = 0. The merger signal starting 
around At = shows excellent agreement for the two configurations and is not affected by 
the spurious signal visible for HD5e at At ~ — SOr^. 

We conclude this discussion with an analysis of the ringdown. After formation of a common 
horizon, the waveform is dominated by an exponentially damped sinusoid, as the merged 
hole rings down into a stationary state. By fitting our results with an exponentially-damped 
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Figure 6.10: / = 2 components of the KI function as generated by a head-on colhsion of black holes 
with initial (coordinate) distance d = 6.37rs (black solid line) and d = lO.STrs (red 
dashed line). The wave functions have been shifted in time such that the formation of 
a common apparent horizon corresponds to At = (and taking into account the time 
it takes for the waves to propagate up to the extraction radius i?ex = BOr^). 



sinusoid, we obtain a characteristic frequency 

rsuj = 0.955 ± 0.005 - i(0.255 ± 0.005) 



(6.4.24) 



This value is in excellent agreement with perturbative calculations, which predict a lowest 
quasinormal frequency r^w = 0.9477 - i0.2561 for / = 2 pMl [T271 [T92] . 

6.4.3.3 Radiated energy 



We now compute the energy flux from the KI master function via equation (5.1.45). The 
fluxes thus obtained for the I = 2 multipole of models HD5e and HD5f in table 6.2 extracted 
at Rex = 60rs, are shown in figure 6.11 As in the case of the KI master function in figure 6.10 



we see no significant variation of the flux for the two different initial separations. The flux 
reaches a maximum value of dE/dt ~ 3.4 x 10~^r5, and is then dominated by the ringdown 
flux. The energy flux from the I = 4 mode is typically four orders of magnitude smaller; this 
is consistent with the factor of 100 difference of the corresponding wave multipoles observed 



in figure 6.9, and the quadratic dependence of the flux on the wave amplitude. Integrating, 
we find that a fraction of E^g^^/M = (8.9 it 0.6) x 10^^ of the centre of mass energy is emitted 
in the form of gravitational radiation. We have verified for these models that the amount of 
energy contained in the spurious radiation is about three orders of magnitude smaller than 
in the physical merger signal. 
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Figure 6.11: Energy flux in the / = 2 component of the KI wave function ^ t, extracted at Rex = 
60rs, for models HD5e (black solid line) and HD5f (red dashed line) in table 6.2 The 
fluxes have been shifted in time by the extraction radius -Rqx = SOrg and the time icAH 
at which the common apparent horizon forms. 



6.5 Discussion 



In this chapter we have presented a framework that allows the generalisation of the present 
generation of 3+1 numerical codes to evolve, with relatively minor modifications, spacetimes 
with SO{D — 2) symmetry in 5 dimensions and SO{D — 3) symmetry in D > 6 dimensions. 
The key idea is a dimensional reduction of the problem, along the lines of references |179|, I184| , 
that recasts the D-dimensional Einstein vacuum equations in the form of the standard four 
dimensional equations plus some source terms. The resulting equations can be transformed 
straightforwardly into the BSSN formulation that has proved remarkably successful in nu- 
merical evolutions of black hole configurations in 3+1 spacetimes. 

The class of problems that may be studied with our framework includes head-on collisions in 
D > 5 and a subset of black hole collisions with impact parameter and spin in D > 6. 

A procedure to construct initial data and a formalism to extract gravitational radiation 
observables from the numerical simulations were also introduced. With these tools, the 
numerical implementation was done by adapting the Lean code and, after a number of 
tests including the convergence of the Hamiltonian and momentum constraints as well as 
comparing numerical results with (semi-) analytic expressions for a single Tangherlini black 
hole in geodesic slicing, we reported results obtained for evolutions of black hole collisions in 
five- dimensional spacetimes. 

As might be expected, stable evolutions of such spacetimes require some modifications of 
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the underlying methods of the so-cahed moving puncture technique, especially with regard 
to the gauge conditions used therein. We have successfully modified the slicing condition 
in order to obtain long-term stable simulations in D = 5 dimensions. Unfortunately, these 
modifications do not appear sufficient to provide long-term stability for arbitrary values of 
the dimensionality D. This issue remains under investigation. 

Besides obtaining the corresponding waveforms for head-on collision of five-dimensional black 
holes, we have further shown that the total energy released in the form of gravitational waves 
is approximately (0.089ib0.006)% of the initial centre of mass energy of the system, for a head- 
on collision of two black holes starting from rest at very large distances. As a comparison, 
the analogous process in D = 4 releases a slightly smaller quantity: (0.055 ± 0.006)%. 

As yet another test of our implementation, the ringdown part of the waveform was also shown 
to yield a quasinormal mode frequency in excellent agreement with predictions from black 
hole perturbation theory. 

The numbers reported here for the total energy loss in gravitational waves should increase 
significantly in high energy collisions, which are the most relevant scenarios for the applica- 
tions described in the Introduction. Indeed, in the four dimensional case, it is known that 
ultra-relativistic head-on collisions of equal mass non-rotating black holes release up to 14% 
of the initial centre of mass energy into gravitational radiation [ID]. The analogous number 
in higher dimensions is as yet unknown, and it remains under investigations using the tools 
here presented. 

Even more energy may be released in high energy collisions with non-vanishing impact 
parameter. In |4H H2] it was shown that this number can be as large as 35% in D = 4. The 
formalism here developed allows, in principle, the study of analogous processes in D > Q. 



6. A Ricci tensor 

In this appendix we give the full details about the computation of the Ricci tensor of 
section 



6.1.1 We start by writing the metric (6.1.1) in block-diagonal form 

ds^ = g^j^e^ = g^,Ax^ dx"" + gij^ (g) , (6.A.1) 

where 

eA = {D^A), 0^ = 8^-6^3% (6.A.2) 
is now the (non-coordinate) basis. Its dual is 

= {dx'', Q'), & = dx' + enB^^dx^". (6.A.3) 

This basis satisfies 



(6.A.4) 
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where J- ab are given by (6.1.8). In the foUowing we will not explicitly assume the expres- 
sions (6.1.8), i.e., we will only assume (unless explicitly mentioned otherwise) that [6^,65] = 

AB^k = J^^ABd%, diQfiv = and ^^Tj — 0- We will not assume the expression for the 
in terms of the coordinate basis (even though we will assume that = 9^). 

Important remark: From now on we will work with the non-coordinate basis ca to simplify 
the calculations. Note that, for an arbitrary tensor Ta, 

Ta\b = deeTA = Cfi {Ta) ^ ObTa = 
In particular, = iD^j = 9^ — CKBjj^di, and as such, 

Tai^ = dejA = (Ta) = d^TA - enE'^diTA. 
We must keep in mind that T4|^^ 7^ 2^a|j^^- On the other hand, eg = 5^ and thus 

= ^kTA. 

We recall our definition of the "covariant derivatives" and Vj as |*] 

Vrpia _ T-) rpia _ , -j^i _r-nloi_ 77/ _'~pioL_ I T~iQf '~n'i^_ t^A rpia,_ (n \ r\ 

VjT-^^ = d-T^"-,^ + r^T^s^ - r'^jT-^, (6.A.6) 
recall also that both connections are metric. 



and that 



We now recall some expressions from section 1.4, which we re- write here for convenience: on 
a non- coordinate basis obeying 

[eA,es] = CAB^dD (6.A.7) 
the torsion-free connection is given by 

r^BC = ]^9^'^ {9db\c + 9dc\b - 9bc\d + CDBC + CDCB - cbcd) , (6.A.8) 
where T^^^c] = —^cbc^, and the Riemann tensor by 

R^BCD = ^'^BD\C ~ ^'^BC\D + Ec'^^ BD — '^^ ED^^ BC — ^^BECCD^ ■ (6. A. 9) 



*Note that now, as we are working on a non-coordinate basis, the order of the indices does matter, i.e., 

r^A. / r",A and r'js / v\j. 
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We use the convention 



Applying these formulas to our case, 

CAB^ = J^AB, cab'" = 0, e^ = D^, ei = di, 



we get from equation (6. A. 8) 







yi _ 


-^^^ 

2'' M*^' 


■pM- — 


2^^ 5jfc-^ Ai/ = ui-i 


"T/^ 

I] 


-25''^Va9jJ' 


-■ = 




f*- = 

m 




f i 


2^*' (s/jife + - 9-fk\l) = r*jfc 



From (6. A. 9) we compute the Ricci tensor 
1 



.A.IO) 



(6.A.11) 



where 



1 



^g^'Vpg-kf^^gij + "-g^'V ^g^-^V^ g.j+ -^g^^gP^gy^g^jT^ p^J"' pa - ^V^V^^ff^j 



(6.A.12) 
(6.A.13) 



1 



,kl\ 



1 



1. 



(6.A.14) 



•r-fc _ _i_ -pfc Tri 

fiu\k ' ^ jk-' 



and we also used T^[ij] = —\c-ij^ = = ^^"^ ~ "S'^M/s" ~ ~2-^"m/3 ~ 0. 

The Ricci scalar is given by 



R = g^'^RAB = g'^'R^u + g^'Rvy 
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We have 

R = R + R- ^g-,jg-^g^-px^TK. " (/^V^g^,-) - li^'g^^V^g^f^ 

- \g'V''^''9-ki^>.gi], (6.A.15) 
where R = g"^^ R-fj- 

We stih need to write the components of the Ricci tensor m the coordinate basis (5^, so 
from here on we need to use the specific form of and, thus, also the algebra in (6.1.8). We 
perform the basis transformation the usual way, 

R = RABe"^ (S) 

Ruu + eKRi^B]^ + eKR^iBi + RjjB'-^Bl^ dx^ (g) dx'' 

Rl^ + eKRijBi) dx* (g) dx'' + (r^j + eKR-fiBf) dx^ (g dx* + ^-dx* (g dx^ 



and thus, in the basis {dfj_,di) where the metric takes the form (6.1.1), we have 

(6.A.16) 



R^i = skR^bI + ^5"^V« {g-i-k^'x,) + \g''ymT9^^^x,m^ + [g^'v^gn 

- (/'V^5s^) = R-^, (6.A.17) 
R^. = R^. + 2eKBlR,)j - e\^Rjr^B^^Bl - \g''^ g-qT\^TK, - Jv, [cpV^^g-i-j 



2 

- |5^^'/'V^5JsV.<7,T - ^V^^^., (6.A.18) 



and 



R = R^R- -^gug'-'^g^^Th.TKu - ,g-,^ - ^'^^ gifJ y^g-Q 

- ^5'V^>'^5fcrV^5^J. (6.A.19) 

Note: We have used the same indices to label components in both the coordinate and non- 
coordinate basis. No confusion shall arise, however, since the non-coordinate basis was 
used merely to simplify the previous computations. In the main text we deal exclusively 
with the coordinate basis, to which all components refer. 



6.B Equations of motion: Einstein frame 



We here write the equations of motion obtained when we write the action (6.1.26) in the 
Einstein frame. To do so, we perform the usual conformal transformation 



giiv = Tp' 



2 

' d-2 



gflu 



(6.B.1) 
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where in our case = e"*^. The action then takes the form 

1 



R + kdf,(l)d^'(j) + n{n - l)e 



-20 



(6.B.2) 



, , d((i-4)(n-l)-2(n+2) , 

where k = „ / -.o and 

{d—2+n)'' 



d-2+n , 



of motion 



y. From action (|6.B.2|) we obtain the equations 
fcV"5„0 + n(n- l)e" 



'2<f> 







n{n-l) 24,- ■ 
d-2 ^ 



(6.B.3) 



Chapter 7 

Non-asymptotically flat spacetimes 



7.1 de Sitter 

Nonlinear dynamics in cosmological backgrounds has the potential to teach us immensely 
about our universe, and also to serve as prototype for nonlinear processes in generic curved 
spacetimes. de Sitter spacetime, as already mentioned in the Introduction, is the simplest 
accelerating universe — a maximally symmetric solution of Einstein's equations with a positive 
cosmological constant — which seems to model quite well the present cosmological accelera- 
tion [Too] . 

Key questions concerning the evolution towards a de Sitter, spatially homogeneous universe 
are how inhomogeneities develop in time and, in particular, if they are washed away by 
the cosmological expansion |193j . Answering them requires controlling the imprint of the 
gravitational interaction between localised objects on the large-scale expansion. Conversely, 
the cosmological dynamics should leave imprints in strong gravitational phenomena like 
primordial black hole formation |194| or the gravitational radiation emitted in a black hole 
binary coalescence, which carry signatures of the cosmological acceleration as it travels 
across the universe. Identifying these signatures is not only of conceptual interest but also 
phenomenologically relevant, in view of the ongoing efforts to directly detect gravitational 
radiation. 

Finally, dynamics in asymptotically de Sitter spacetimes could also teach us about more 
fundamental questions such as cosmic censorship: two black holes of sufficiently large mass 
in de Sitter spacetime would, upon merger, give rise to too large a black hole to fit in its 
cosmological horizon. In this case the end state would be a naked singularity. This possibility 
begs for a time evolution of such a configuration. Does the time evolution of non-singular 
data containing two black holes result in a naked singularity, or are potentially offending 
black holes simply driven away from each other by the cosmological expansion? 

In this section, following [ 195j . we report on numerical evolutions of black hole binaries 
in an asymptotically de Sitter geometry. Even though we consider a range of values for the 
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cosmological constant far larger than those which are phenomenologicahy viable, these results 
provide useful insight on the general features of dynamical black hole processes in spacetimes 
with a cosmological constant, which can improve our understanding of our universe. 



7.1.1 Evolution equations 



The Einstein equations with cosmological constant A are 

1 



R 



:gnuR — —J^Qiiu , 



(7.1.1) 



and we will always consider A > 0. We perform the 3+1 decomposition by introducing 
the projection operator 7^,^ and the normal to the three dimensional hyper-surface S, n'^ 
{n^Uf^ = —1), as outlined in section 2.5 and write the evolution equations in the BSSN 



form (4.1.14) 



From (7.1.1), we straightforwardly compute the source terms 



Svr^; = A , 



8ttS = -3A . 



(7.1.2) 



A new evolution variable x = exp(2iyA/3t)x has been introduced instead of the usual BSSN 
variable x |193j . The reason is that for black hole evolutions it is crucial to impose a floor 
value on x, typically 10~^ or 10~^, which is inconsistent with the natural behaviour of this 
variable in a de Sitter spacetime (as we will see below): x~^ ~ exp(2y^A/3t). In contrast 
X — >■ 1 when r —>■ 00 for all times. The evolution equations are thus 



dtx = 2x{aK - di^/^ + P'dix + m 3 X : 



dtK 
dtAij 



aA . 



(7.1.3a) 

(7.1.3b) 

(7.1.3c) 
(7.1.3d) 
(7.1.3e) 



where [•••] denotes the right-hand side of the BSSN equations (4.1.14) in the absence of 
source terms. 



7.1.2 Schwarzschild-de Sitter 



The Schwarzschild-de Sitter spacetime, solution of (7.1.1), written in static coordinates reads 

-2 - -f{R)dT^ + f{R)-^dR^ + R^dn2 . (7.1.4) 



The solution is characterised by two parameters: the black hole mass m and the Hubble 
parameter H, 

(7.1.5) 



f{R) = 1 - 2m/ R - H^R^ , H = y^A/3 
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f{R) has two zeros, at i? = R±, < if 

0<mH < mH^rit , mH^rit = Vl/27 . (7.1.6) 

These zeros are the location of the black hole event horizon (R-) and of a cosmological 
horizon (i?+). U H = 0, then R_ = 2m; if m = 0, then R+ = 1/H. U H,m 0, then 
R^ > 2m and i?+ < 1/H. Since R is the areal radius, the area of the spatial sections of the 
cosmological horizon decreases in the presence of a black hole; and the area of the spatial 
sections of the black hole horizon increases in the presence of a cosmological constant, as one 
would intuitively anticipate. 

The basic dynamics in this spacetime may be inferred by looking at radial timelike geodesies. 
They obey the equation (di?/dr)^ = E'^ — f{R), where r is the proper time and E is the 
conserved quantity associated to the Killing vector field d/dT. In the static patch (i?_ < 
R < i?+), E can be regarded as energy. From this equation we see that f{R) is an effective 
potential. This potential has a maximum at 

i?ma. = (m/i/2) 1/3. (7.I.7) 

Geodesies starting from rest (i.e. di?/dT(r = tq) = 0) will fall into the black hole if R- < 
R < i?max or move away from the black hole if i?max < R < R+- 

As we will discuss in the next section, the initial data for an evolution in the de Sitter universe 
can be computed in a similar manner as has been done in asymptotically flat space as long 
as one chooses a foliation with extrinsic curvature Kij having only a trace part. Such a 
coordinate system is known for Schwarzschild-de Sitter: McVittie coordinates |l96j . These 
are obtained from static coordinates by the transformation (T, R) — )• (t, r) given by 

r 7? d /? 

R = {l + ^fa{t)r , T = t + H ^ ^ ^ , (7.1.8) 



fiR)y^l-2m/R 

where a{t) = exp{Ht) and ^ = 2c!Jt)r • obtains McVittie's form for Schwarzschild-de 

Sitter: 



ds' = -(^^^J dt' + a{ty{l + (f{dr' + rMn2) . (7.1.9) 
For t = constant, one can show that indeed Kj = —H5j. 

By setting m = in McVittie coordinates one recovers an FRW cosmological model with 
k = (flat spatial curvature) and an exponentially growing scale factor. The cosmological 
horizon Tic discussed above, located at R = 1/H, stands at r-^^ = l/{He^*). The spatial 
sections of 1-Lc seem to be shrinking down in this coordinate system. What happens, in 
fact, is that the exponentially fast expansion is taking any observer to the outside of T-Lc- 
This is a well known phenomenon in studies of inflation and, as we shall see, has important 
consequences for the numerical evolution. 
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7.1.3 Numerical Setup 

The cosmological constant introduces a new term (when compared with the vacuum case) in 



the Hamiltonian constraint obtained after the canonical 3+1 decomposition (2.6.3c) 



R - KijK'^ + K'^ = 2A (7.1.10) 



In references \197\ I198j it was observed that imposing a spacetime shcing obeying 
—HSj, and a spatial metric of the form dP = ijj^jij dx'^dx^ , the equations to be solved in 
order to obtain initial data are equivalent to those in vacuum. In particular, for a system of 
black holes momentarily at rest (with respect to the given spatial coordinate patch), the 
conformal factor ip takes the form 

N 

^ = i + y (7.1.11) 

There are + 1 asymptotically de Sitter regions, as |r — r(j)| — t- 0, +oo; the total mass for 
observers in the common asymptotic region (|r — r^j^l — t- +oo) is [X98\. 

Boundary conditions for all quantities are imposed by looking at the behaviour of massless 
perturbing fields in a pure de Sitter background. Accordingly, we impose the following 
asymptotic behaviour for all BSSN variables 

dtf - dtfo + -^drf + ^—^ -Hif-fo) = 0. (7.1.12) 
a[tj a[t)r 

We should note that we also performed evolutions using different sets of boundary conditions, 
to test the independence of the results on boundary conditions imposed in a region with no 
causal contact with the interaction region. As far as the behaviour and location of the 
horizons and all quantities discussed in this paper are concerned, no noticeable difference 
could be found. 



Our numerical simulations use the Lean code |149| . see section 4.3 The calculation of Black 
hole Apparent Horizons (BAHs) and Cosmological Apparent Horizons (CAHs) is performed 
with AHFinderDirect |153| 1154] . We remark that BAHs, found as marginally trapped 
surfaces, indicate in de Sitter space (with the same legitimacy as in asymptotically flat space) 
the existence of an event horizon |199j . CAHs are surfaces of zero expansion for ingoing null 
geodesies. In a single black hole case, in McVittie coordinates, the black hole event horizon 
and cosmological horizon are indeed foliated by apparent horizons. 

The "expanding" behaviour of the coordinate system led us to add a new innermost refinement 
level at periodic time intervals so as to keep the number of points inside the cosmological 
horizon approximately unchanged. The necessity for adding extra refinement levels effectively 
limits our ability to follow the evolution on very long timescales, as the number of time steps 
to cover a fixed portion of physical time grows exponentially. This feature resembles in many 
ways the recently reported work by Pretorius and Lehner on the follow-up of the black string 
instability 
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7.1.4 Numerical Results 

As a first test of tlie numerical implementation, we performed evolutions of a single black 



hole imposing the McVittie slicing condition; that is, we use (7.1.9) as initial data and impose 
dta = AmrHe^*/{m + 2re^* f , 



0. 



(7.1.13) 



throughout the evolution. The analytical solution (7.1.9) can be compared with the numerical 
results. For a single black hole evolution with m = 1 and H = 0.8-ffcrit) the results are 
displayed in figure 7.1 Using this slicing, the runs eventually crash (at t ~ 12m). By 
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=0 


■ 
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= 5m 


• 
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= 10m 




Figure 7.1: Conformal factor x for a single black hole evolution with H = O.Sffcrit using the McVittie 

results are plotted, along 
0), against the expected 



slicing condition, equation (7.1.13). The obtained numerical results are plotted, along 

x{z) imposed at z 



the z coordinate (symmetry xi^z) 
analytical solutions (solid lines). 



contrast, the standard "1+log" slicing condition (4.2.3) 



dta = (3'dia -2a{K - Kq) 



(7.1.14) 



where Kq = —3H = — v3A, enables us to have long term stable evolutions. As consistency 
checks, the areal radii at the apparent horizons (both black hole horizon and cosmological 
horizon) are constants in time and have the value expected from the analytical solution in 
a single black hole spacetime. Moreover, the areal radius at fixed coordinate radius evolves 
with time in the way expected from the exact solution. 



For binary black hole initial data, we start by reproducing the results of Nakao et. al |198| . 
where the critical distance between two black holes for the existence of a common BAH 
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already at t = was studied. We thus prepare initial data (7.1.11) with mi 
all quantities in units of the total mass m 



m2 and take 



nil + "^2 ■ The two punctures are set initially 
at symmetric positions along the z axis. The critical value for the cosmological constant, for 
which the black hole and cosmological horizon coincide is now mHcnt = 1/ \/27- We call small 
(large) mass binaries those, for which H < -fTcrit (H > -fTcrit)- Our results for the critical 
separation in small mass binaries, at t = 0, as function of the Hubble parameter are shown in 
figure 7.2 The line (diamond symbols) agrees, after a necessary normalisation, with figure 14 
of [T98] . 




Figure 7.2: Critical coordinate distance for small mass binaries, from both initial data and dynamical 
evolutions, as well as a point particle estimate, as a function of H/Hcru- We obtain this 
estimate from the coordinate distance to the horizon, equation (7.1.7), for a particular 
value of m. The t — line refers to the critical separation between having or not having 
a common BAH in the initial data. The inset shows details of the approach to the 
critical line for H = 0.6-ffciit, where a is an acceleration parameter. 



We now consider head-on collisions of two black holes with no initial momentum, i.e. the time 
evolution of these data. We have monitored the Hamiltonian constraint violation level for 
cases with and without cosmological constant. We observe that the constraint violations are 



comparable in the two cases and plot in figure 7.3 a snapshot of the Hamiltonian constraint 



violation at t = 48m for parameters H = 0.9-ffcrit and d = 0.8m, a typical case with non-zero 
cosmological constant. We have used two resolutions, m/160 and m/192 (on the innermost 
refinement level) and have rescaled the dashed curve by Q2 = (192/160)2 as expected for 
second-order convergence. 

For subcritical Hubble constant H < Hcnt = l/(\/27"i), we monitor the evolution of the 
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z/m 

Figure 7.3: Hamiltonian constraint violation along the z-axis at time t = 48m for a simulation with 
H — 0.9i?crit and initial distance d = 0.8m. 



areal radius of the BAHs and that of the CAH of an observer at z = 0. For instance, for 
H = 0.9-ffcrit and proper (initial) separation 3.69?7i we find that the areal radii of the BAH 
and CAH are approximately constant and equal to Rbau ~ 2.36m and -RcAH — 4.16m, 
respectively. As expected the two initial BAHs, as well as the final horizon, are inside the 
CAH. As a comparison, a Schwarzschild-de Sitter spacetime with the same H has Rbah — 
2.43m and i?cAH — 4.16m. This suggests that the interaction effects (binding energy and 
emission of gravitational radiation) are of the order of a few per cent for this configuration. 

As the initial separation grows, so does the total time for merger. For separations larger than 
a critical value, the two black holes do not merge, but scatter to infinity. For such scattering 
configurations, the simulations eventually exhibit a regime of exponentially increasing proper 
distance between the BAH. Just as in scatters of high energy black holes [l2], here we find 
that the immediate merger/scatter regimes are separated by a blurred region, where the holes 
sit at an almost fixed proper distance for some time; cf. figure [774} By performing a large set 
of simulations for various cosmological parameters H and initial distance d, we have bracketed 
the critical distance for the merger /scatter region as a function of the Hubble parameter H 
for the "dynamical" case, i.e., the initial coordinate distance between the black holes such 
that no common BAH forms. The results are displayed in figure 7.2 (circles and x symbols). 



As expected the critical distance becomes larger as compared to the initial data value ( "t = 0" 
line): there are configurations for which a common BAH is absent in the initial data but 
appears during the evolution (just as in asymptotically flat spacetime). The numerical results 
can be qualitatively well approximated by a point particle prediction — from equation (7.1.7). 
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Figure 7.4: Proper distance between the black hole horizons as a function of time for the H = 
0.9iJcrit, and initial (coordinate) distance d ~ 0.9m. The two holes stay at approximately 
constant distance up to t ^ 8m after which cosmological expansion starts dominating. 



To do such comparison a transformation to McVittie coordinates needs to be done; we have 
performed such transformation at McVittie time t = 0. Intriguingly, for a particular value 
of m ~ 0.7, the point particle approximation matches quantitatively very well the numerical 
result; the curve obtained from the geodesic prediction in figure [7^ is barely distinguishable 
from the numerical results. 

A further interesting feature concerns the approach to the critical line. For an initially static 
binary close to the critical initial separation, the coordinate distance d scales as d = do + at'^. 
In general the acceleration parameter scales as log a = C + Tlog{d — do), where T = 1 in 
the geodesic approximation. A fit to our numerical results for H = 0.6-ffcrit (dashed curve in 



the inset of figure 7.2) for example yields C = —3.1, T = 0.9 in rough agreement with this 



expectation. Details of this regime are given in the inset of figure 7.2 



Finally, we have performed evolutions with H > Hcnt- On the assumption of weak grav- 
itational wave release, such evolutions can test the cosmic censorship conjecture since the 
observation of a merger in such case would reveal a violation of the conjecture |200j . From 
general arguments and from the simulations with H < i/criti we know the cosmological 
repulsion will dominate for sufficiently large initial distance and in that case we can even 
expect that a CAH for the observer at z = will not encompass the BAHs. This indicates 
the black holes are no longer in causal contact and therefore can never merge. Our numerical 
results confirm this overall picture. To test the potentially dangerous configurations, we focus 
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on the regime in which the black holes are initially very close. A typical example is depicted 
in figure [T^l for a supercritical cosmological constant H = 1.05-ffcrit; and an initial coordinate 
distance d/m = 1.5002. Even though the initial separation is very small, we find that the 
holes move away from each other, with a proper separation increasing as the simulation 
progresses. In fact, further into the evolution, a distorted CAH appears, and remains for as 
long as the simulation lasts. At late times, this CAH is spherically symmetric, and has an 
areal radius which agrees, to within 10~^, with that of an empty de Sitter spacetime with the 
same cosmological constant. The evolution therefore indicates that the spacetime becomes, 
to an excellent approximation, empty de Sitter space for the observer at z = and that the 
black holes are not in causal contact. Observe that qualitatively similar evolutions can be 
found in small mass binaries when the initial distance is larger than the critical value 




-0.6 0.0 0.6 -0.6 0.0 0.6 -0.6 0.0 0.6 



z/m z/m z/m 

Figure 7.5: Snapshots at different times (from left to right t/m — 0.0, 8.0156, 20.016) of a simulation 
with H = 1.05i/crit: and an initial coordinate distance d/m = 1.5002. The dotted blue 
line denotes the CAH (for an observer at z = 0) which is first seen in this simulation 
at t/m — 8.0156, highly distorted. At late times, the CAH has an areal radius of 
R = 4.94876, while the "theoretical value" for pure dS is i? = 1/7? = 4.94872, a 
remarkable agreement showing that the spacetime is accurately empty dS for the observer 
at z = and the black holes are not in causal contact. 



7.1.5 Final Remarks 

We have presented evidence that the numerical evolution of black hole spacetimes in de 
Sitter universes is under control. Our results open the door to new studies of strong field 
gravity in cosmologically interesting scenarios. In closing, we would like to mention that 
our results are compatible with cosmic censorship in cosmological backgrounds. However, 
an analytic solution with multiple (charged and extremal) black holes in asymptotically de 
Sitter spacetime is known, and has been used to study cosmic censorship violations |201j . In 
collapsing universes a potential violation of the conjecture has been reported, although the 
conclusion relied on singular initial data. To clarify this issue, it would be of great interest to 
perform numerical evolution of large mass black hole binaries, analogous to those performed 
herein, but in collapsing universes. This will require adaptations of our setup, since the 
"expanding" behaviour discussed of the coordinate system will turn into a "collapsing" one. 
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which raises new numerical challenges. 



7.2 Black holes in a box 

Anti-de Sitter (AdS) is a non-globally hyperbolic spacetime, which essentially means that 
it is not enough to prescribe a set of evolution equations and some initial configuration in 
order to predict what will happen in the future. On such spacetimes, the boundary plays 
an active role, and in order to have a well-defined Cauchy problem the initial data (and 
evolution equations) must be supplemented by appropriate boundary conditions at the time- 
like conformal boundary. 

In this section, we will give a brief overview of the work presented in |104j where a "toy 
model" for Anti-de Sitter was considered by imprisoning a black hole binary in a box with 
mirror-like boundary conditions and thus exploring the active role that boundary conditions 
play in the evolution of a bulk black hole system. 



7.2.1 Numerical setup 



The vacuum Einstein equation are written in the BSSN scheme (4.1.14), introduced in 
section 4.1, with the gauge conditions (4.2.3), (4.2.6). We evolve these equations with the 



Lean code [149J, see section 4.3 



The work here presented differs from previous implementations of the Lean code (and most 
other codes) in the treatment of the outer boundary conditions, which is herein considered 
to be a refiecting sphere or, rather, an approximation of it by using so-called Lego spheres; 
cf. section 3 in |202j . The numerical implementation of such boundary conditions is schemat- 



ically illustrated in figure 7.6 



Points outside the outer circle of radius Rb + ^R are not required for updating regular points 
and are simply ignored in the numerical evolution. In practice, we ensure that the boundary 
shell is always of sufficient thickness to accommodate discretization stencils required for the 
update of regular gridpoints. The specific boundary condition is then determined by the 
manner in which we update grid functions on the boundary points marked as x in the figure. 

To mimic the global structure of an AdS spacetime we thus enclose the black hole binary 
inside the spherical mirror and set 

|/ = 0, (7.2.1) 



at each boundary point with / denoting any of the BSSN variables listed in equations (4.1.14) 



The final ingredient needed for our numerical implementation is related with the spurious 
radiation present when evolving black hole binaries — often called junk radiation — which can 
be traced back to the methods used to compute the initial data. To avoid contamination of 
our simulations by such spurious radiation being trapped inside our reflective boundary we 
employ standard outgoing radiation boundary conditions at early times and only switch on 
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Figure 7.6: Illustration of a (Lego-)spherical outer boundary. 



our reflective condition at 

*ref = T^B + Atpulse, (7.2.2) 

where we estimate the duration of the spurious wave pulse Atpuisc from previous simulations 
of similar setups in asymptotically flat spacetimes as for example presented in |203| 11491 1204| . 
The spurious radiation is thus given sufficient time to leave the computational domain. 

Wave extraction is employed in the fashion outlined in section |5.1.1[ where we will herein 
also measure the ^'o Weyl scalar, which encodes the incoming gravitational wave signal. 

The results we will report in the following sections all refer to an inspiral simulation with 
total mass M = M\ +M2, where the black hole punctures were set with an initial coordinate 
distance of d = 6.517M and with Bowen-York momentum parameter Pj = ib0.133M. The 
grid structure used was 

{(48, 24, 12, 6) X (1.5, 0.75), h = 1/56} , (7.2.3) 
and the Weyl scalars have been extracted at r^x = 35M. 
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7.2.2 Gravitational wave signal 



The nature of our specific configuration is ideal to study both the outgoing (^'4) as well as 
the ingoing (^0) gravitational wave pulses. 



The gravitational wave signal is dominated by the quadrupole contributions which is shown 
in figure 7.7 The ingoing signal xp22 been shifted in time by At = 10 M to compensate 
for the additional propagation time from the extraction radius rex = 35 M to the boundary 
Rb = 4:0 M and back after reflection. The reflection introduces an additional phase shift of 
Ai;^ = TT which has also been taken into account in the figure. Within numerical errors, we 
find the resulting outgoing and subsequent ingoing pulses to overlap. 
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Figure 7.7: Real part of the I = m = 2 mode of rM'^Q and rM'^4. The ingoing signal rMvfo has 
been shifted in time by At — lOM and in phase by tt (thus equivalent to an extra minus 
sign) to account for the additional propagation time and the reflection. 



7.2.3 Interaction of the wave pulse with the remnant black hole 

We define black hole mass in terms of the equatorial radius of the horizon Cg by |205) 

M = ^ . (7.2.4) 

Ait 



In figure 7.8 it is shown the fractional deviation (M — Mq)/Mq of the mass of the final black 
hole from its value immediately after merger together with the irreducible mass and the black 
hole spin J . The mass remains approximately constant until the pulse returns after its first 
refiection, then increases, remains constant during the second passage of the pulse and so on. 
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In contrast, the spin only shows a significant increase during the first scattering of the pulse 
off the black hole. 




Comparing the increase in the horizon mass with the amount of gravitational wave energy 
radiated during the last stages of the inspiral, plunge and merger of a corresponding binary 
system in an asymptotically flat spacetime — which is about 3.5% of the total energy of the 
system |204l I149j — we estimate that about 15% of the energy emitted during the merger is 
absorbed by the central spinning black hole per interaction. 



7.2.4 Final remarks 

In this section, we have given just a brief overview of the work presented in |lU4j where the 
global structure of an AdS background was mimicked by introducing a reflecting wall at a 
finite radius. Inside this cavity a black hole inspiral was evolved. 

The results presented are consistent with the intuitive expectations for a wavepacket of 
radiation (generated during inspiral plus merger) travelling back and forth between the 
mirror-like wall and the black hole: part of this radiation is absorbed when interacting with 
the black hole (especially high-frequencies). We estimate that about 15% of the wavepacket's 
energy is absorbed by the black hole per interaction, at least during the first cycles. 

It would be extremely interesting to extend this work to implement the evolution of black 
holes in real AdS backgrounds, following the recent works in |6 H 162^ [63]. 



CHAPTER 7. NON-ASYMPTOTICALLY FLAT SPACETIMES 



110 



7.3 Black holes in cylinders 

From the gauge/gravity duality to braneworld scenarios, black holes in compactified space- 
times play an important role in fundamental physics. Our current understanding of black 
hole solutions and their dynamics in such spacetimes is rather poor because analytical tools 
are capable of handling a limited class of idealised scenarios, only. 

In this section, following |206j . we wish to study how the compactness of extra dimensions 
changes the dynamics of such higher-dimensional gravity scenarios. There is considerable 
literature on Kaluza-Klein black holes and black holes on cylinders |106| I107| I105| I108| : 
the full non-linear dynamics of black holes on such spacetimes, however, seems to remain 
unexplored. 



7.3.1 Setup 

We are interested in describing the evolution of black holes in a five dimensional spacetime 
with one periodic direction. For a five dimensional cylindrical Minkowski spacetime, M^'^ x , 
the metric can be written as 



-di^ + dx^ + dy^ + 



+ dz^ 



(7.3.1) 



51 



The direction is parameterised by z, which takes values in the interval [— L,L], with the 
two endpoints identified and L G M"*". The coordinate also parameterises a circle, this circle 
is, however, homotopic to a point, since it shrinks down to zero size at y = 0, where y is a 
radial coordinate in the y — (j) plane which is part of M^'^ — figure 7.9 




Figure 7.9: Illustration of the coordinate system for the Minkowski spacetime M^''^ x S^. A slice 
with t = constant and x = constant is shown, y, (p parameterise a plane, wherein y is 
a radial direction and an azimuthal coordinate. At each point in this plane there is 
a non-contractible circle parameterised by z. This is illustrated by exhibiting this circle 
on various points along an orbit of d/d(j) and also at y = 0. Space-time is a (trivial) 
bundle over M^^'^. 



Following the approach outlined in section 6.2 we take our five dimensional metric ansatz to 
be 

ds^ = g^.udx^dx" + \{x^')d(t)'^ , (7.3.2) 
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where x^^ = {t,x,y,z). We perform a dimensional reduction by isometry on and end up 
with a four dimensional model of gravity coupled to a scalar field. Performing the standard 
3 + 1 decomposition and writing the equations in the BSSN scheme, the evolution equations 
of the resulting system are those of (4.1.14) with matter terms given by (6.2.23). We here 
use periodic boundary conditions along the z direction and Sommerfeld radiative boundary 
conditions along x and y. 



7.3.2 Initial data 



Following the approach of section 6.3, the four-dimensional Brill-Lindquist initial data ap- 
propriate to describe non-spinning, non-rotating black holes momentarily at rest, take the 
form 



-fijdx'dx^ = ^P'^[dx^ + dy^ + dz^] , A = 



Ki, =0 = Kx. 



In a spacetime with standard topology (wherein z parameterises a line), the initial data for 
two black holes with horizon radius r^^ and punctures placed at {x,y,z) = (0,0,iba), takes 
the form 



+ 



4[x2 +y'^ + {z - a)2] 4[x2 + y'^ + {z + a^] 



(7.3.3) 
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'z = -L 



r 



-2L 



Single black hole 



t 



z = 2L + a 
z — 2L — a 

~z = L 

^z — a 
iz — —a 

'z = -L 

z = -2L + 1 
z = -2L - < 
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Figure 7.10: Illustration of the correspondence between real space and covering space for a single 
black hole (left panel) and a situation of head on collision (right panel). The dashed 
boxes drawn in the covering space contain a single copy of the real space setup and 
correspond also to what is contained in the numerical grid. 

The appropriate initial data to describe a black hole in can be viewed as having an infinite 
array of black holes, all with the same mass, separated by coordinate distance Az = 2L — 



figure 7.10 Since the superposition of various black holes in a line is described by adding up 
the corresponding initial data, for the infinite array of two black holes in the circle located 
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at z = ±a (0 < a < L) with horizon radii r^, i = 1,2 (or, equivalently, for two black holes in 
S^) the initial data is given by 

^ 4[x2 Jry'^ + {z-a- 2Ln)^] ^ 4[x2 + y'i + {z + a - 2Ln)2] 

(7.3.4) 



?1= — oo 

^ , vr(r^)^ sinhf ^ vr(r|)2 sinh f 



SLp cosh f - cos 8Lp cosh ^ - cos ' 

where = + and in the last equality we have used the result in |105j . 



7.3.3 Results 



Again, we use the Lean code, introduced in section 4.3, for the numerical evolutions. 
The main difference to standard implementations of Lean is the use of periodic boundary 
conditions, which is also non-trivial to implement in a parallel code. 

We will now show some results obtained for a head-on collision (from rest) of black holes with 
an initial separation of 10.37 rs, i.e., a = 5.185 rg- The z {z ^ [— L, L]) coordinate has been 
compactified with L/rs = 64, 32, 16. For comparison purposes, we have also performed 
a simulation with "standard" outgoing boundary conditions (L — )• oo), which will be here 
referred to as "outgoing" . 

All results will be presented in units of the Schwarzschild radius rg = rs,i + rs,2- 



7.3.3.1 Hamiltonian constraint 



Figure 7.11 shows the Hamiltonian constraint along the x and z axis, respectively, for several 
time steps for the L/rs = 32 case. As we can see, the constraint is indeed being satisfied 
with high accuracy. 





20 40 



100 120 



20 25 



Figure 7.11: Hamiltonian constraint along the x and z axis, for the L/rs — 32 case. 
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7.3.3.2 Collision time 



Next, we study the changes in cohision time for different compactification radii. Whereas 
we do not observe any (noticable) difference for L/rs = 64 and L/rs = 32 as compared to 
the outgoing case, the case L/rs = 16 shows already a noticable difference. The puncture 



trajectories for these cases are plotted in figure 7.12 



Puncture trajectory 




Figure 7.12: Puncture z coordinate as function of time for the outgoing and L/rg — 16 cases. 



Recall that one can think of a black hole in a cylindrical space as an array of infinite black 
holes. Therefore, for a head-on collision on such a cylindrical space, each black hole will also 
feel the gravitational pull of all the other black holes. Naively, one thus expects that for this 
cylindrical case it will take longer for the black holes to collide, which is what we observe in 



figure 7.12 



7.3.4 Final remarks 



Using the formalism introduced in section 6.1 , we were able to reduce the head-on collision of 
(non-spinning) black holes on cylindrical spacetimes (in any dimension) to an effective 3 + 1 
system with a scalar field, and used this procedure to successfully evolve a head-on collision 
of two black holes on a five-dimensional cylindrical spacetime. 

Further issues that we wish to investigate include monitoring the deformation of the black 



holes' apparent horizon and computing the energy radiated, along the lines of section 6.4.2.3 



We also further plan to perform simulations with smaller compactification radii and study 
the equivalent six-dimensional system. 



Chapter 8 



Einstein-Maxwell 



In this last chapter we go back to four dimensions once again, this time in Einstein-Maxwell 
theory, to perform fully non-linear numerical simulations of charged black hole collisions |207| . 

As mentioned in the Introduction, the dynamics of binary systems of charged, i.e. Reissner- 
Nordstrom (RN), black holes have remained rather unexplored territory. Perhaps this is due 
to the expectation that astrophysical black holes carry zero or very small charge; in particular, 
black holes with mass M, charge Q and angular momentum aM^ are expected to discharge 
very quickly if Q/M > W-^^{a/M)-^/'^{M/MQy/^ [2081 [Il2]. There is, nevertheless, a good 
deal of motivation for detailed investigations of the dynamics of charged black holes. 

In the context of astrophysics, charged black holes may actually be of interest in realistic 
systems. First, a rotating black hole in an external magnetic field will accrete charged particles 
up to a given value, Q = 2BqJ Thus it is conceivable that astrophysical black holes 

could have some (albeit rather small) amount of electrical charge. Then it is of interest 
to understand the role of this charge in the Blandford-Znajek mechanism |112j . which has 
been suggested for extracting spin energy from the hole, or in a related mechanism capable 
of extracting energy from a moving black hole |110^ I113j to power outflows from accretion 
disk- fed black holes. Numerical simulations of charged black holes interacting with matter 
and surrounding plasma will enable us to study such effects. 

Motivation for the numerical modelling of charged black holes also arises in the context 
of high energy collisions. It is expected that trans-Planckian particle collisions form black 
holes; moreover, well above the fundamental Planck scale such processes should be well 
described by general relativity and other interactions should become negligible [(PO, an idea 
poetically stated as matter does not matter for ultra high energy collisions [66j. But is this 
expectation really correct? Calculations of shock wave collisions suggest that even though 
other interactions — say charge — may become irrelevant in the ultra- relativistic limit, the 
properties of the final black hole (and of the associated emission of gravitational radiation) 
do depend on the amount of charge carried by the colliding particles |114t I115j . This issue 
can be clarified by the simulation of high-energy collisions of charged black holes and the 
subsequent comparison of the results to those obtained for electrically neutral systems. 



114 
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Finally we note a variety of conceptual aspects that merit a more detailed investigation of 
charged black hole systems. In head-on collisions with small velocity, the intuition borrowed 
from Larmor's formula in Minkowski space suggests a steady growth of the emitted power with 
the acceleration. However, it is by now well established that for uncharged black holes the 
gravitational radiation strongly peaks near around time of formation of a common apparent 
horizon. Does the electromagnetic radiation emission follow a similar pattern? And what is 
the relative fraction of electromagnetic to gravitational wave emissions? Moreover, a non- 
head on collision of charged non-spinning black holes will allow us to study, as the end 
state, a (perturbed) Kerr-Newman geometry, which would be extremely interesting: linearised 
perturbations around Kerr-Newman black holes do not decouple |209^ 1192) and so far close 
to nothing is known about their properties. Among others, the stability of the Kerr-Newman 
metric is an outstanding open issue. Furthermore, it has been observed that the inspiral 
phase of an orbiting black-hole-binary system can be well understood via post-Newtonian 
methods |210j (see also e.g. [231 I211j ). The additional radiative channel opened by the 
presence of electric charge provides additional scope to probe this observation. 

With the above motivations in mind we here initiate the numerical study of non-linear 
dynamics of binary systems of charged black holes, building on previous numerical evolutions 
of the Einstein-Maxwell system |116l 11091 \in\ I118j . For reasons of simplicity, we focus in 
this study on binary systems for which initial data can be constructed by purely analytic 
means |133l I212j : head-on collisions, starting from rest, of non-spinning black holes with 
equal charge-to-mass ratio. This implies in particular that the black holes carry a charge 
of the same sign, so the electromagnetic force will always be repulsive. We extract both 
gravitational and electromagnetic radiation and monitor their behaviour as the charge-to- 
mass-ratio parameter of the system is varied. 



8.1 Evolution equations 

We will adopt the approach outlined in [213^ I117j to evolve the electro-vacuum Einstein- 
Maxwell equations which incorporates suitably added additional fields to ensure the evolution 
will preserve the constraints. This amounts to considering an enlarged system of the form 



where •kF^'^ denotes the Hodge dual of the Maxwell- Faraday tensor F^'^, k is a constant 
and the four-velocity of the Eulerian observer. We recover the standard Einstein-Maxwell 
system of equations when ^ = = With the scalar field ^ and pseudo-scalar $ introduced 
in this way, the natural evolution of this system drives ^ and $ to zero (for positive k), thus 
ensuring the magnetic and electric constraints are controlled \213\ I116j . The electromagnetic 




8ttT, 



(8.1.1) 
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stress-energy tensor takes the usual form 

1 

Ait 



.1.2) 



8.1.1 3+1 decomposition 



We employ the 3 + 1 decomposition, as explained in section 2.5 where we introduced the 
3-metric 

l^J.v = Qixv + n^j.nu , (8.1.3) 

and further decompose the Maxwell- Faraday tensor into the more familiar electric and mag- 
netic fields measured by the Eulerian observer moving with four velocity 

where we use the convention €1230 = \/—9, ^a(5-y = ^aP'ySn^ ■, ^123 = ^/l- 



.1.4) 



We write the evolution equations in the BSSN form (4.1.14) where, for the case of the 
electromagnetic energy-momentum tensor of equation (8.1.2), the source terms are given 

by 

E = T^'^n^n, = {E'Ei + B'Bi) , 



1 

47r 



(8.1.5) 



1 

47r 



-EiEj - BiBj + -jij ( E^'Ek + B'^Bk 



and S = Y'' Sij- The evolution of the electromagnetic fields is determined by equation (8.1.1 ) 
whose 3+1 decomposition becomes 



{dt - Cp) E' = aKE' + e'^'^x-' [%iB'd,a + a [B^d.^i + ^udjB' - x'^lkiB'djX 

- axf^dj^ , 

{dt - Cp) B' = aKB' - e'^\-^ [ikiE'd^a + a (E^dj^i + ^udjE' - x'^iE'Ojx 

- axf^dj^ , 

{dt - Cp) ^ 
{dt - £/?) $ 



-aViE^ - aK^ , 



(8.1.6) 

Here, denotes the Lie derivative along the shift vector /3*. The Hamiltonian and momen- 
tum constraint are 

n = R + K'^ - K'^Kij - WttE = , 

3 2 (8.1.7) 

Mi = DjA^ - -A^x~'djX - ^diK - 8vrj- = , 

where Di is the covariant derivative associated with the three- metric jij. 
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8.1.2 Initial data 

We focus here on black hole binaries with equal charge and mass colliding from rest. For 
these configurations, it is possible to construct initial using the Brill-Lindquist construction 



(outlined in section 3.2.1 for the vacuum spacetimes; see |133tl212] for the charged case). The 



main ingredients of this procedure are as follows. 

For a vanishing shift time symmetry implies Kij = 0. Combined with the condition of 
an initially vanishing magnetic field, the magnetic constraint DiB^ = and momentum con- 
straint are automatically satisfied. By further assuming the spatial metric to be conformally 
flat 

-fijdx'dx^ = (dx^ + dy2 + dz"^) , (8.1.8) 
the Hamiltonian constraint reduces to 

Ai/j = -^E^iP^ , (8.1.9) 

where A is the flat space Laplace operator. The electric constraint. Gauss's law, has the 
usual form 

DiE' = . (8.1.10) 

Quite remarkably, for systems of black holes with equal charge-to-mass ratio, these equations 
have known analytical solutions |212j . For the special case of two black holes momentarily at 
rest with "bare masses" mi, m2 and "bare charges" qi, q2 = qim2/rai this analytic solution 
is given by 

V 2|x-xi| 2\x-X2\J A\\x-xi\ \x-x2\y (8 111) 

E —ip [qi -z ^-pr + q2 -p ^-fo- 

\ \X — Xi\-^ \X — X2\'^ 

where Xi is the coordinate location of the ith. "puncture" Q 

The initial data is thus completely specified in terms of the independent mass and charge 
parameters mi, m2, qi and the initial coordinate separation d of the holes. These uniquely 
determine the remaining charge parameter q2 via the condition of equal charge-to-mass ratio. 
In this study we always choose mi = m2 and, without loss of generality, position the two holes 
symetrically around the origin such that zi = d/2 = —Z2- The resulting initial three metric 



7ij follows from equations. (8.1.8), (8.1.11) while the extrinsic curvature Kij and magnetic 



field vanish on the initial slice. 

We use the same gauge conditions and outer boundary conditions for the BSSN variables as 



used in vacuum simulations, cf. equations (4.2.3) and (4.2.6). As outer boundary condition 
for the electric and magnetic fields we have imposed a falloff as 1/r^ — from (8.1.11). For 



the additional scalar fields a satisfactory behaviour is observed by imposing a falloff as 1/r^ 
(which is the expected falloff rate from dimensional grounds). 



*We note that this fohation, in isotropic coordinates, only covers the outside of the external horizon. 
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8.2 Wave Extraction 



For a given set of initial parameters mi = m2, Qi = q2, d, the time evolution provides us 
with the spatial metric 'jij, the extrinsic curvature Kij as well as the electric and magnetic 
fields as functions of time. These fields enable us to extract the gravitational and 

electromagnetic radiation as explained in section 5.1.1 Details concerning the numerical 
implementation can be found in [149] . 



We recall that the radiated flux and energy are given by the expressions (5.1.21 ) and (5.1.28): 



Feu 



dt 

d^EM 

dt 



lim 



[ dtv""(t'] 

l.m -^-"^ 



lim — 

r— >oo 4:71 



Lm 



(8.2.1) 
(8.2.2) 



As is well known from simulations of uncharged black-hole binaries, initial data obtained 
from the Brill-Lindquist construction contains "spurious" radiation, which is an artifact of 
the conformal-flatness assumption. In calculating properties of the radiation, we account 



for this effect by starting the integration of the radiated flux in equations (8.2.1), (8.2.2) at 



some finite time At after the start of the simulation, thus allowing the spurious pulse to first 
radiate off the computational domain. In practice, we obtain satisfactory results by choosing 
At = Rex + 50 M. Because the physical radiation is very weak for both the gravitational 
and electromagnetic channel in this early infall stage, the error incurred by this truncation 



is negligible compared with the uncertainties due to discretization; cf. section 8.4.4 



8.3 Analytic predictions 

Before discussing in detail the results of our numerical simulations, it is instructive to discuss 
the behaviour of the binary system as expected from an analytic approximation. Such an 
analysis not only serves an intuitive understanding of the binary's dynamics, but also provides 
predictions to compare with the numerical results presented below. 

For this purpose we consider the electrodynamics of a system of two equal point charges in 
a Minkowski background spacetime. As in the black hole case, we denote hy qi = q2 = Q/2 
and mi = m2 = M/2 the electric charge and mass of the particles which are initially at rest 
at position z = ibd/2. 

It turns out useful to first consider point charges in Minkowski spacetime in the static limit. 
The expected behaviour of the radial component of the resulting electric field is given by |214j 



(8.3.1) 
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which for a system of two charges of equal magnitude at z = ibd/2 becomes 



Ef ~ 




(8.3.2) 



The dipole vanishes in this case due to the reflection symmetry across z = 0. This symmetry 
is naturally preserved during the time evolution of the two-charge system. Furthermore, the 
total electric charge Q is conserved so that the leading-order behaviour of the electromag- 
netic radiation is given by variation of the electric quadrupole, just as for the gravitational 
radiation. Notice that in principle other radiative contributions can arise from the accel- 
erated motion of the charged black holes. From experience with gravitational radiation 
generated in the collision of electrically neutral black-hole binaries, however, we expect this 
"Bremsstrahlung" to be small in comparison with the merger signal and hence ignore its 
contributions in this simple approximation. The good agreement with the numerical results 
presented in the next section bears out the validity of this quadrupole approximation. In 
consequence, it appears legitimate to regard the "strength" of the collision and the excitation 
of the black-hole ringdown to be purely kinematic effects. 

An estimate for the monopole and quadrupole amplitudes in the limit of two static point 
charges is then obtained from inserting the radial component of the electric field (8.3.2) into 



the expression (5.1.25) for <I>i and its multipolar decomposition (5.1.26) 



rVf = \l "—Qd^ ~ 0.59Qfi^ . 



TiQ « 1.71Q, 




(8.3.3) 
(8.3.4) 



The expectation is that these expressions provide a good approximation for the wave sig- 
nal during the early infall stage when the black holes are moving with small velocities. 

Ji after the merger and 
should eventually approach zero as a single merged 



Equation (8.3.3) should also provide a good approximation for 



ringdown whereas the quadrupole 
hole corresponds to the case d = in equation (|8.3.4|). 



In order to obtain analytic estimates for the collision time and the emitted radiation, we 
need to describe the dynamic behaviour of the two point charges. Our starting point for this 
discussion is the combined gravitational and electromagnetic potential energy for two charges 
z = 1, 2 in Minkowski spacetime with mass and charge m-j, qi at distance r from each other 



V 



Gmim2 



+ 



1 qiq2 



(8.3.5) 



r 47reo r 

For the case of two charges with equal mass and charge rrii = M/2, q^ = Q/2 and starting 
from rest at zq = ib(i/2, conservation of energy implies 



Mz^ 



where we have used units with G = 47ren 



M'^B _ 
Az 

= 1 and 



M'^B 
2d 



(8.3.6) 



B = l- Q^jM'^ 



(8.3.7) 
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The resulting equation of motion for z{t) is obtained by differentiating equation (8.3.6) which 
results in 

Mz = -|^ + ^ = -M2A. (8.3.8) 

8z^ 8z^ 82;^ 

An estimate for the time for collision follows from integrating equation (8.3.6 ) over z £ [d/2, 0] 

/ ^collision \ TT d 



(8.3.C 



Prom the dynamic evolution of the system we can derive an approximate prediction for 
the electromagnetic radiation by evaluating the (traceless) electric quadrupole tensor Qij = 
J d?xp{x){2>XiXj — r'^Sij) [214j . In terms of this quadrupole tensor, the total power radiated 
is given by |214j 

ft« = Ei^5(ib<3«- (8-3.10) 

For clarity we have reinstated the factors 47reo and c' here. Using 

-^{z'^) = 6iz + 2z'z , (8.3.11) 



and the equations of motion (8.3.6), (8.3.8) we find 

Fem . (8.3.12) 

Using J dt{. . .) = J dz/z{. . .), we can evaluate the time integral up to some cutoff separation, 
say Zmin = Cihb, where b is the horizon radius of the initial black hole, b = M(l + \fB)l2 and 
afe = 0(1) is a constant. This gives, 

^ = B'^'Hl^l^C^ ~ 2a,6)3/2(15d2 + 2Ma,b + ?,2ay) 
M ^ 50400((ia66)7/2 ' ^ ■ ■ ) 

Emission of gravitational radiation follows from the quadrupole formula, which is a numerical 
factor 4 times larger, and where the charge is be replaced by the mass, 

^ = r5/2^7/2 {d-2a,bfl\lhd'' + 2Ma,b + 2,2alb^) 
M 12600((iQb6)7/2 ■ y ■ ■ ) 

For Q = 0, CKfe = 1, d = 00 we thus obtain 

= — 0.0012 , (8.3.15) 

M 840 ' ^ ' 



in agreement to within a factor of 2 with numerical simulations (see |176j and table 8.1 below; 
the agreement could be improved by assuming Of, ~ 1.3). As a general result of this analysis 
we find in this approximation. 

For non-extremal holes Q < M, our analytic considerations therefore predict that the energy 
emitted in electromagnetic radiation is at most 25% of the energy lost in gravitational 
radiation. As we shall see below, this turns out to be a remarkably good prediction for 
the results obtained from fully numerical simulations. 



CHAPTER 8. EINSTEIN-MAXWELL 



121 



8.4 Numerical Results 



The numerical integration of the Einstein-Maxwell equations (4.1.14), (8.1.6) has been per- 
formed using fourth-order spatial discretisation with the Lean code, originally presented 



in |149j for vacuum spacetimes, see section 4.3 



The initial parameters as well as the grid setup and the radiated gravitational and electromag- 
netic wave energy for our set of binary configurations is listed in table 8.1 All binaries start 
from rest with a coordinate distance d/M ~ 8 or d/M ~ 16 while the charge-to-mass ratio 
has been varied from Q/M = to Q/M = 0.98. Note that identical coordinate separations 
of the punctures for different values of the charge Q/M correspond to different horizon-to- 
horizon proper distances. This difference is expected and in fact analysis of the RN solution 
predicts a divergence of the proper distance in the limit Q/M — )• 1. 



Table 8.1: Grid structure in the notation of section II E of [1491. coordinate distance c?/M, proper 
horizon-to-horizon distance L/M , charge Q/M, gravitational [E^^) and electromagnetic 
i-^rad ) radiated energy for our set of simulations. The radiated energy has been computed 
using only the / = 2, to = mode; the energy contained in higher-order multipoles such 
as Z = 4, TO = is negligible for all configurations. 



Run 


Grid 




d/M 


L/M Q/M 


^rad 








dOSqOO 


{(256,128,64,32,16,8) 


X (2, 1,0.5), 1/80} 


8.002 


11.56 





5.1 X 10" 


-4 






d08q03 


{(256,128,64,32,16,8) 


X (2, 1,0.5), 1/80} 


8.002 


11.60 


0.3 


4.5 X 10" 


-4 


1.3 X 10" 


5 


d08q04 


{(256,128,64,32,16,8) 


X (2, 1,0.5), 1/80} 


8.002 


11.65 


0.4 


4.0 X 10" 


-4 


2.1 X 10" 


5 


d08q05c 


{(256,128,64,32,16,8) 


X (2, 1,0.5), 1/64} 


8.002 


11.67 


0.5 


3.3 X 10" 


-4 


2.7 X 10" 


5 


d08q05m 


{(256,128,64,32,16,8) 


X (2, 1,0.5), 1/80} 


8.002 


11.70 


0.5 


3.4 X 10" 


^4 


2.7 X 10" 


5 


d08q05f 


{(256,128,64,32,16,8) 


X (2, 1,0.5), 1/96} 


8.002 


11.67 


0.5 


3.4 X 10" 


^4 


2.7 X 10" 


5 


d08q055 


{(256,128,64,32,16,8) 


X (2, 1,0.5), 1/80} 


8.002 


11.70 


0.55 


3.0 X 10" 


-4 


2.89 X 10 


-5 


d08q06 


{(256,128,64,32,16,8) 


X (2, 1,0.5), 1/80} 


8.002 


11.75 


0.6 


2.6 X 10" 


-4 


2.97 X 10 


-5 


d08q07 


{(256,128,64,32,16,8) 


X (2, 1,0.5), 1/80} 


8.002 


11.87 


0.7 


1.8 X 10" 


^4 


2.7 X 10" 


5 


d08q08 


{(256,128,64,32,16,8) 


X (2, 1,0.5), 1/80} 


8.002 


12.0 


0.8 


9.8 X 10" 


-5 


1.8 X 10" 


5 


d08q09 


{(256,128,64,32,16,8) 


X (2, 1,0.5), 1/80} 


8.002 


12.3 


0.9 


2.6 X 10" 


-5 


5.5 X 10" 


6 


d08q098cc 


{(256,128,64,32,16,8) 


X (2, 1,0.5), 1/64} 


8.002 


12.3 


0.98 


7.0 X 10" 


-7 


2.1 X 10" 


7 


d08q098c 


{(256,128,64,32,16,8) 


X (2, 1,0.5), 1/80} 


8.002 


13.1 


0.98 


4.3 X 10" 


-7 


1.4 X 10" 


7 


d08q098m 


{(256,128,64,32,16,8) 


X (2, 1,0.5), 1/96} 


8.002 


13.1 


0.98 


3.4 X 10" 


-7 


1.0 X 10" 


7 


d08q098f 


{(256,128,64, 32,16,8) 


X (2, 1,0.5), 1/112} 


8.002 


13.0 


0.98 


4.0 X 10" 


-7 


9.5 X 10" 


8 


d08q098fr 


{(256,128,64, 32,16,8) 


X (2, 1,0.5), 1/128} 


8.002 


13.0 


0.98 


4.05 X 10 


-7 


8.75 X 10 


-8 


d08q098fff {(256,128,64,32,16,8) 


X (2, 1,0.5), 1/136} 


8.002 


13.1 


0.98 


3.73 X 10 


-7 


8.41 X 10 


-8 


dl6q00 


{(256,128,64,32,16) x 


(4, 2, 1,0.5), 1/64} 


16.002 


20.2 





5.5 X 10" 


-4 






dl6q05 


{(256,128,64,32,16) x 


(4, 2, 1,0.5), 1/64} 


16.002 


20.3 


0.5 


3.6 X 10" 


-4 


2.9 X 10" 


5 


dl6q08 


{(256,128,64,32,16) x 


(4, 2, 1,0.5), 1/80} 
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-4 
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5 


dl6q09 


{(256,128,64,32,16) x 


(4, 2, 1,0.5), 1/80} 


16.002 


21.0 


0.9 


2.7 X 10" 


-5 


5.9 X 10" 
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8.4.1 Code tests 



Before discussing the obtained results in more detail, we present two tests to validate the 
performance of our numerical implementation of the evolution equations: (i) single black- 
hole evolutions in geodesic slicing which is known to result in numerical instabilities after 
relatively short times but facilitates direct comparison with a semi-analytic solution and (ii) 
convergence analysis of the radiated quadrupole waveforms for simulation d08q05 of table [STT} 



The geodesic slicing condition is enforced by setting the gauge functions to a = 1, /3* = 
throughout the evolution. The space part of the Reissner-Nordstrom solution in isotropic 



coordinates is given by equation (8.1.8) with a conformal factor [215\ I2l6] 



1 + 



M 
2r 



47-2 



^.4.11 



The time evolution of this solution is not known in closed analytic form, but the resulting 
metric components can be constructed straightforwardly via a simple integration procedure. 
As expected, we find a time evolution in this gauge to become numerically unstable at times 
r of a few M. Before the breaking down of the evolution, however, we can safely compare 
the numerical and "analytical" solutions. This comparison is shown in figure |8.1| for the ^zz 
component of the spatial metric and the component of the electric field and demonstrates 
excellent agreement between the semi-analytic and numerical results. 
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Figure 8.1: The numerical profiles for and (symbols) obtained in geodesic slicing at various 
times r are compared with the semi-analytic results (lines). 



For the second test, we have evolved model d08q05 using three different resolutions as 
listed in table 8.1 and extracted the gravitational and electromagnetic quadrupole (/ = 
2, m = 0) at i?ex = lOOM. For fourth-order convergence, we expect the differences between 
the higher resolution simulations to be a factor 2.78 smaller than their coarser resolution 
counterparts. The numerically obtained differences are displayed with the corresponding 
rescaling in figure |8.2[ Throughout the physically relevant part of the waveform, we observe 
the expected fourth-order convergence. Only the spurious initial radiation (cf. the discussion 
at the end of section 8.2) at early times At < —20 in the figure exhibits convergence closer to 
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second order, presumably a consequence of high-frequency noise contained in this spurious 
part of the signal. From Richardson extrapolation of our results we estimate the truncation 
error of the radiated waves to be about 1%. The error due to extraction at finite radius, on 
the other hand, is estimated to be 2 % at Rex = lOOM. 




20 40 60 

At/M 



3i 







— '^c-K 






-- 2.78(h„, -hj) . 







20 40 60 80 100 120 

At/M 



Figure 8.2: Convergence analysis for simulation d08q05 of table 8.1 with resolutions he — Af/64, 
hjn — M /SO and hf = M/96. The panels show differences of the (2, 0) multipoles of the 
real parts of 5*4 (left) and $2 (right) extracted at i?ex = 100 M; in each case, the high- 
resolution differences have been rescaled by a factor 2.78 as expected for fourth-order 
convergence. 



8.4.2 Collisions of two black holes: the "static" components and infall time 



We start the discussion of our results with the behaviour of the gravitational and electro- 
magnetic multipoles when the system is in a nearly static configuration, i.e. shortly after the 
start of the simulation and at late stages after the ringdown of the post-merger hole. At 



these times, we expect our analytic predictions (8.3.3), (8.3.4) for the monopole and dipole 



of the electromagnetic field to provide a rather accurate description. Furthermore, the total 
spacetime charge Q is conserved throughout the evolution, so that the monopole component 



of $1 should be described by (8.3.3) at all times. The quadrupole, on the other hand, is 



expected to deviate significantly from the static prediction (8.3.4) when the black holes start 
moving fast. 



As demonstrated in figure 8.3 we find our results to be consistent with this picture. Here 
we plot the monopole and quadrupole of $1. The monopole part (left panel) captures 
the Coulomb field and can thus be compared with the total charge of the system. It is 
constant throughout the evolution to within numerical error and shows agreement with the 



analytic prediction of equation (8.3.3) within numerical uncertainties; we measure a slightly 



smaller value for the monopole field than expected from the total charge of the system, but 
the measured value should increase with extraction radii and agree with the total charge 
expectation at infinity. This is consistent with the extrapolation of the measured value to 
infinity as shown in the figure. The quadrupole part (right panel) starts at a non-zero value in 
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Figure 8.3: Monopole 0^ (left) and quadrupole (right) of the radial part of the elect roraagnetic 
field $1 extracted at i?cx = lOOA/ for simulation d08q05 of table |8.l] The dashed curves 

oo in the static limit. For the 



show the predictions of equations (8.3.3), (8.3.4) at R 



monopole case, we also added the curves obtained by extrapolating the results to infinite 
extraction radius; these curves — dotted lines — essentially overlap with the predictions 
from equation 



.3.3) 



excellent agreement with equation (8.3.4), deviates substantially during the highly dynamic 
plunge and merger stage and eventually rings down towards the static limit cpf^ = as 
expected for a spherically symmetric charge distribution. 



The analytic approximation of section 8.3 also predicts a value for the time of collision (8.3.9) 



for a given set of initial parameters. In particular, we see from this prediction that for fixed 
initial separation d and mass M the collision time scales with the charge as icoiUsion ~ 1/ V^- 
In comparing these predictions with our numerical results we face the difficulty of not having 
an unambiguous definition of the separation of the black holes in the fully general relativistic 
case. From the entries in table [8?T] we see that the proper distance L varies only mildly for 
fixed coordinate distance d up to Q/M « 0.8. For nearly extremal values of Q, however, L 
starts increasing significantly as expected from our discussion at the start of this section. We 
therefore expect the collision time of the numerical simulations rescaled by VS/to, where to 
is the corresponding time for the uncharged case, to be close to unity over a wide range of 
Q/M and show some deviation close to Q/M = 1. This expectation is borne out in figure 8.4 



where we show this rescaled collision time, determined numerically as the first appearance of 
a common apparent horizon, function of Q/M. 



8.4.3 Waveforms: infall, merger and ringdown 



The dynamical behaviour of all our simulations is qualitatively well represented by the 



waveforms shown in figure 8.5 for simulations dlGqOO, dl6q05 and dl6q09. The panels show 
the real part of the gravitational (left) and electromagnetic (right) quadrupole extracted at 
-Rex = 100 M as a function of time with At = defined as the time of the global maximum of 



the waveform. From the classical analysis (8.3.10), we expect the waveforms \I'4, $2 to scale 



roughly with B and the mass or charge of the black holes (the scaling with B is non-trivial. 
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Figure 8.4: Time for apparent horizon formation, re-scaled by the factor \/B and the apparent 
horizon formation time to for an electrically neutral binary. We note that the change in 
the quantity we plot is only, at most, of 2%. The coordinate time itself, however, varies 
by a factor 5 as one goes from Q = to Q — 0.98M. 



but both an analytic estimate and the numerical results indicate the scaling is approximately 
linear, which we shall therefore use for re-scaling the plots in the figure). 
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Figure 8.5: Real part of the (2, 0) mode of 5*4 (left) and $2 (right panel) extracted at Rex — lOOM. 



The early stage of the signals are marked by the spurious radiation due to the construction 
of initial data which we ignore in our analysis. Following a relatively weak phase of wave 
emission during the infall of the holes, the radiation increases strongly during the black- 
hole merger around At = in the figure and decays exponentially as the final hole rings 
down into a stationary state. This overall structure of the signals is rather similar for 
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the electromagnetic and the gravitational part and follows the main pattern observed for 
gravitational- wave emission in head-on collisions of uncharged black holes |176t 1217] , 

Table 8.2: Comparison of the ringdown frequencies obtained from (i) perturbative calculations |192j 
and (ii) fitting a two-mode profile to the numerically extracted waveforms. For Q/M = 
the electromagnetic modes are not excited. For values of Q/AI > 0.9 the electromagnetic 
mode becomes so weak that we can no longer unambiguously identify it in the numerical 
data. 
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- 0.096i 
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0.382 - 


0.0896f 


0.381 


- 0.09H 




0.494 - 
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The final, exponentially damped ringdown phase is well described by perturbation tech- 
niques |192j . In particular, charged black holes are expected to oscillate with two different 
types of modes, one of gravitational and one of electromagnetic origin. For the case of 
vanishing charge, the electromagnetic modes are not present, but they generally couple 
for charged black holes, and we expect both modes to be present in the spectra of our 
gravitational and electromagnetic waveforms. For verification we have fitted the late-stages 
of the waveforms to a two-mode, exponentially damped sinusoid waveform 

f{t) = Aie-*"i* + A2e-''^^\ (8.4.2) 

where Ai are real-valued amplitudes and coi complex frequencies. The results are summarised 



in table 8.2 for selected values of the charge-to-mass ratio of the post-merger black hole. Real 
and imaginary part of the fitted frequencies agree within a few percent or better with the 
perturbative predictions. For the large value Q/M, however, the wave signal is very weak 
and in such good agreement with a single ringdown mode (the gravitational one) that we 
cannot clearly identify a second, electromagnetic component. This feature is explained once 
we understand how the total radiated energy is distributed between the gravitational and the 
electromagnetic channels. For this purpose, we plot in figure [STBl the Fourier spectrum of the 
relevant wavefunctions or, more precisely, their dominant quadrupole contributions obtained 
for simulation d08q03 IV'^^'P, where for any function / 

/oo 
e'^^'/m. (8.4.3) 
-oo 

It is clear from the figure that most of the energy is carried in the fundamental gravitational- 
wave like mode with a peak at approximately u ~ 0.37, close to the oscillation frequency of 



the fundamental gravitational ringdown mode; see table 8.2 
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Figure 8.6: Power spectrum for the gravitational (long dashed) and electromagnetic (short dashed) 
quadrupole extracted from simulation d08q03. Note that the spectrum peaks near the 



fundamental ringdown frequency of the gravitational mode; cf. table 8.2 



8.4.4 Radiated energy and fluxes 



The electromagnetic and gravitational wave fluxes are given by equations (8.2.1) and (8.2.2). 



We have already noticed from the waveforms in figure 8.5 that the electromagnetic signal 



follows a pattern quite similar to the gravitational one. The same holds for the energy flux 
which is shown in figure 8.7 for a subset of our simulations with Q/M = 0, 0.5 and 0.9. 
From the figure, as well as the numbers in table |8.1[ we observe that the energy carried by 
gravitational radiation decreases with increasing Q/M, as the acceleration becomes smaller 



and quadrupole emission is suppressed, in agreement with prediction (8.3.14) 



This is further illustrated in figure |8.8[ which illustrates the radiated energy carried in 
the gravitational quadrupole and the electromagnetic quadrupole as well as their ratio as 
functions of the charge-to-mass ratio Q/M. For the case of vanishing charge, the total 
radiated energy is already known from the literature; e.g. jl76] . The value increases mildly 
with the initial separation as a consequence of the slightly larger collision velocity but is 
generally found to be close to E^^ /M = 0.055%. Our values of 0.051% for d/M ~ 8 and 
0.055% for d/M ~ 16 are in good agreement with the literature. As we increase Q/M, 
however, E^^ decreases significantly and for Q/M = 0.9 (0.98) has dropped by a factor 
of about 20 (10^) relative to the uncharged case. For practical reasons, we have explored 
the largest ratio Q/M = 0.98 for the smaller initial separation d/M ~ 8 only; the near 
cancellation of the gravitational and electromagnetic interaction and the resulting slow-down 
of the collision lead to a very long infall stage with essentially zero dynamics. 



In contrast to the monotonically decreasing gravitational-wave energy, the electromagnetic 
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Figure 8.7: Radiated fluxes for simulations d08q05, d08q09 and dOSqOO of table 8.1 We have aligned 
the curves in time such that their global maximum coincides with t = 0. The inset shows 
the exact same plot with the y-axis in logarithmic units. 
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Figure 8.8: Energy radiated in the gravitational and electromagnetic quadrupole as well as the ratio 
of the two as a function of Q/M. 
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signal reaches a local maximum around Q/M = 0.6, an expected observation as the elec- 
tromagnetic radiation necessarily vanishes for Q/M = (no charge) and Q/M = 1 (no 
acceleration) but takes on non-zero values in the regime in between. Closer analysis of our 



classical, flat-space calculation (8.3.13) predicts a maximum electromagnetic radiation output 
at 



.4.4) 



Qmax = \l " M « 0.605M , 

in excellent agreement with the results of our simulations. 

We finally consider the ratio of electromagnetic to gravitational wave energy (dotted curve 
in figure |8.8|) 




As predicted by our analytic calculation 



.3.16), this ratio increases mono- 
tonically with Q/M for fixed separation d. A fit of our numerical results yields E^^/E^^ = 
0.27 Q^/M'^ and for our largest value Q/M = 0.98, we obtain a ratio of 0.227 to be compared 
with ~ 0.24 as predicted by equation (8.3.16). Bearing in mind the simplicity of our analytic 
model in section f 



3l the quantitative agreement is remarkable. 



8.5 Conclusions 

In this chapter, we performed a numerical study of collisions of charged black holes with 
equal mass and charge in the framework of the fully non- linear Einstein-Maxwell equations. 
Our first observation is that the numerical relativity techniques (formulation of the evolution 
equations, gauge conditions and initial data construction) developed for electrically neutral 
black hole binaries can be straightforwardly extended to successfully model charged binaries 
even for nearly extremal charge-to- mass ratios Q/M < 1. In particular, we notice the 
contrast with the case of rotating black holes with nearly extremal spin which represents 
a more delicate task for state-of-the-art numerical relativity; cf. references [27\ [28] for the 
latest developments on this front. This absence of difficulties for charged holes is not entirely 
unexpected. Considering the construction of initial data, for instance, an important difference 
arises in the customary choice of conformally flat Bowen-York initial data [134] which greatly 
simplifles the initial data problem. While the Kerr solution for a single rotating black 
hole does not admit conformally flat slices |136| and therefore inevitably results in spurious 
radiation, especially for large spin parameters, this difficulty does not arise for charged, but 



non-rotating black holes; cf. equation (8.4.1) and |215| . 



The excellent agreement between the classical calculation for the energy emission and the 
numerical results reported here, allow for an investigation of cosmic censorship close to 
extremality. If we take two black holes with Mi = M2 = M/2, Q^ = Q2 = (M - 5)/2 
and we let them fall from infinity, to first order in 6 we get 

Qtot = M -6 

^ . (8.5.1) 

Mtot = M - ^rad 



Now, the classical result (8.3.14) implies that the dominant term for the radiated energy is 
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£^rad ~ B^/^M ~ Thus we get 




^ 1 - ^ 




(8.5.2) 



where A; is a constant. We conclude that cosmic censorship is preserved for charged colhsions 
of nearly extremal holes {6 <^ M), on account of the much longer collision time, which yields 
much lower velocities and therefore much lower energy output. The differences between the 
cases of spinning mergers and charged collisions are interesting. In the former case, naked 
singularities are avoided by radiation carrying away more angular momentum (via orbital 
hangup [34])- In the latter case, our results suggest that naked singularities are avoided by 
the smaller radiation emission, due to the smaller accelerations involved in the infall. 

We have here evolved a sequence of binaries, with equal charge-to-mass ratio starting from 
rest, with Q/M varying from zero to values close to extremality. Starting with the electrically 
neutral case, where our gravitational wave emission E^^ /M = 0.055% agrees well with the 
literature, we observe a monotonic decrease of the emitted gravitational wave energy as 
we increase Q/M. For our largest value Q/M = 0.98, E'^^ is reduced by about three 
orders of magnitude, as the near cancellation of the gravitational and electromagnetic forces 
substantially slows down the collision. In contrast, the radiated electromagnetic energy 
reaches a maximum near Q/M = 0.6 but always remains significantly below its gravitational 
counterpart. Indeed, the ratio Ef^/E^^ increases monotonically with Q/M and approaches 
about 25% in the limit Q/M — )• 1. We find all these results to be in remarkably good qual- 
itative and quantitative agreement with analytic approximations obtained in the framework 
of the dynamics of two point charges in a Minkowski background. This approximation also 
predicts that the collision time relative to that of the uncharged case scales ~ ^1-QVM2 
which is confirmed within a few percent by our numerical simulations. 



Chapter 9 



Final remarks 



This probably just goes to show something, 
but I sure don't know what. 

Calvin 
Calvin & Hobbes 

Numerical relativity is a fantastic tool to study and explore spacetimes whose exact form is 
not known. 

After decades of efforts, the first stable, long-term evolutions of the orbit and merger of two 
black holes were finally accomplished in 2005, and since then considerable progress has been 
made. This field has now reached a state of maturity, and several codes and tools exist that 
allow one to perform evolutions of black holes — with quite generic initial configurations — in 
standard four dimensional vacuum gravity. 

In addition to the original (main) motivation coming from the two-body problem, it was 
quickly realised that numerical relativity could be helpful for a much broader range of 
scenarios, with some motivation coming from fields other than gravity itself. 

In this work we have thus worked to extend numerical relativity tools to new frontiers, 
opening a range of uncharted territory in black hole physics to be explored with contemporary 
numerical relativity. In particular, we have presented the following: 

(i) a dimensional reduction procedure that allows the use of existing 3 + 1 numerical 
codes to evolve higher-dimensional spacetimes with enough symmetry, including head- 
on collisions in Z) > 5 and black hole collisions with impact parameter and spin in 
D>6; 

(ii) a generalisation of the TwoPunctures spectral solver, allowing for the computation of 
initial data for a boosted head-on collision of black hole binaries in higher-dimensional 
spacetimes; 

(iii) a wave extraction procedure that allows the extraction of gravitational radiation observ- 
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ables from numerical evolutions of head-on collisions of black holes in D dimensions; 

(iv) with the above tools, numerical simulations of black hole collisions from rest in five- 
dimensional spacetimes were successfully evolved, the corresponding wave forms were 
obtained and total energy released in the form of gravitational waves was computed; 

(v) evolutions of black holes in non-asymptotically flat spacetimes, including asymptotically 
dc Sitter spacetimes, "boxed" spacetimes with mirror-like boundary conditions, and 
five-dimensional cylindrical spacetimes; 

(vi) numerical evolutions of collisions of charged black holes with equal mass and charge, and 
a calculation of the energy released via emission of gravitational and electromagnetic 
radiation. 

Several open questions and research avenues remain to be explored, and we thus close with 
a list of natural sequels for this program: 

• A systematic investigation of black hole collisions and dynamics in generic dimension. 
Even though the formalism here presented is valid in arbitrary dimension, the long-term 
numerical stability of the implementation is a different matter altogether. Currently, 
only the five-dimensional case seems to be relatively robust, with numerical instabil- 
ities occurring in all D > 5 cases tried so far. It is possible that such instabilities 
may be cured with a suitable choice of gauge conditions. These issues remain under 
investigation. 

• Related to the previous point, it could be of interest to systematically investigate the 
merits and disadvantages, from the point of view of the numerical implementation, 
of dimensional reduction procedures (such as the one here presented) versus evolution 
schemes that make use of the Cartoon method. 

• The numbers here reported for the total energy loss for the five-dimensional black hole 
head-on collisions refer to collisions from rest. For the applications described in the 
Introduction, however, high velocity collisions are the most relevant ones. Such cases 
do not seem to be as robust as the analogous four dimensional systems, with numerical 
instabilities appearing when large boost parameters are considered. Investigation on 
this front is still under way. 

• For the Einstein-Maxwell study, a natural step is considering more generic types of 
initial data, in order to tackle some of the issues discussed in the Introduction. A 
non-zero boost, for instance, will allow us to study both binary black hole systems 
that will coalesce into a Kerr-Newman black hole and the impact of electric charge 
on the dynamics and wave emission (electromagnetic and gravitational) in high energy 
collisions. A further interesting extension is the case of oppositely charged black holes. 
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